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SUMMARY 


Southwest  Research  Institute  conducted  a  research  project  focused  on  developing  a  methodology 
for  calculating  the  probabilistic  response  of  deep  underground  tunnels  as  part  of  an  ongoing 
program  sponsored  by  DNA.  A  specific  objective  of  the  effort  was  to  develop  a  methodology 
that  is  more  efficient  than  traditional  Monte  Carlo.  Prior  to  this  project,  probabilistic  analysis  of 
underground  structures  was  made  using  Monte  Carlo  simulation;  however,  due  to  it's  inherent 
inefficiency,  only  analytical  or  very  simple  numerical  tunnel  models  could  be  employed. 
Recently,  it  has  been  well  established  by  the  DNA  geomechanics  community  that  predicting 
tunnel  response  to  highly  dynamic  loading  requires  highly-detailed  numerical  models,  which  in 
turn  require  substantial  computational  resources  and  time. 

The  major  result  of  this  project  is  the  development  and  verification  of  a  suite  of  probabilistic 
analysis  methods  that  allow  all  significant  uncertainties  in  complex  numerical  tunnel 
deformation/damage  models  to  be  simulated  in  an  efficient  manner.  Specifically,  new  and 
enhanced  existing  probabilistic  analysis  algorithms  have  been  developed  and  integrated  with  a 
general-purpose  transient  dynamic  finite  element  program.  We  participated  in  two  large 
verification  and  validation  programs  moderated  by  DNA  to  gain  confidence  in  the  underlying 
deterministic  tunnel  models.  Several  demonstration  probabilistic  analyses  were  performed  and 
verified  using  Monte  Carlo  and/or  other  simulation  methods.  In  all  cases,  the  advanced 
probabilistic  methods  used  were  shown  to  be  as  accurate  as  Monte  Carlo  simulation,  but  orders 
of  magnitude  more  efficient. 

Several  significant  advances  have  been  made  within  the  probabilistic  algorithms  to  tailor  their 
use  for  underground  tunnel  analysis:  a  capability  for  computing  efficiently  the  confidence 
bounds  on  the  calculated  probability  of  failure  due  to  random  and  systematic  errors;  a  capability 
for  handling  non-normal  correlated  random  variables;  new  methods  for  performing  sensitivity 
analysis  to  assess  the  impact  of  changing  the  mean,  standard  deviation,  or  distribution  type  of 
any  input  random  variable;  and  a  methodology  for  treating  the  estimated  coefficients  in  the  finite 
element  constitutive  model  as  random  variables. 

Specific  accomplishments  include: 

1 .  An  analytical  model  for  predicting  tunnel  closure  was  formulated, 
implemented,  and  integrated  with  SwRI's  fast  probability  integration 
(NESSUS/FPI)  program.  Probabilistic  sensitivity  factors  computed  as  a 
by-product  of  the  FPI  calculation  provided  new  insight  into  the 
importance  ranking  of  the  input  variables. 

2.  A  general  automated  interface  was  developed  between  the  probabilistic 
algorithms  in  NESSUS/FPI  and  the  PRONTO  explicit  dynamic  finite 
element  (FE)  program.  Any  variable  pertaining  to  the  analysis  can  be 
considered  random,  e.g.,  material  parameters,  loading  time  history,  and 
boundary  conditions.  Probabilistic  methods  such  as  advanced  mean  value, 
adaptive  importance  sampling,  standard  Monte  Carlo,  first  order  reliability 
method,  second  order  reliability  method,  and  fast  probability  integration 
can  all  be  used. 
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3.  A  new  method  for  estimating  confidence  bounds  on  the  tunnel  closure 
cumulative  distribution  function  was  developed  and  implemented  in  FPL 
Confidence  bounds  are  used  to  estimate  the  uncertainty  band  in  a 
probabilistic  prediction  due  to  effects  such  as  statistical  uncertainty,  bias 
and  systematic  errors.  The  new  method  is  based  on  FPI's  advanced  mean 
value  (AMV)  probabilistic  algorithm.  Before  this  method  was  developed, 
Monte  Carlo  simulation  was  required  to  estimate  the  confidence  bounds. 

4.  A  new  and  general  method  for  treating  non-normal  correlated  random 
variables  in  the  probabilistic  analysis  was  developed  and  implemented  in 
FPL  This  capability  allows  the  correlation  coefficient  between  any  of  the 
random  variables  to  be  input  regardless  of  the  distribution  type  specified 
for  each  random  variable.  Prior  to  this  development,  only  random 
variables  having  certain  distribution  functions  could  be  correlated. 

5.  A  probabilistic  cap  model  methodology  was  developed  and  implemented. 

The  development  of  the  model  focused  on  the  accurate  description  of  the 
parameter  uncertainties  in  the  constitutive  model.  The  current  model 
treats  the  parameters  as  random  variables,  since  they  can  be  measured  and 
statistically  modelled.  Statistical  parameters  are  estimated  from  repeat 
laboratory  tests. 

6.  Several  sensitivity  calculations  have  been  developed  and  implemented. 
First-order  approximate  sensitivities  with  respect  to  the  mean  and  standard 
deviation  of  each  input  random  variable  are  computed  from  the 
probabilistic  sensitivity  factors  produced  as  a  by-product  of  the  FPI 
calculation.  A  new  method  for  computing  these  same  sensitivities  which 
re-uses  the  samples  used  by  the  adaptive  importance  sampling  (AIS) 
method  has  been  developed.  Finally,  the  sensitivity  of  the  calculated 
probability  with  respect  to  the  choice  of  input  distribution  has  also  been 
developed. 

7.  Several  demonstrative  probabilistic  analyses  have  been  performed  using 
1)  an  analytical  tunnel  closure  model;  2)  a  transient  dynamic  FE  model  of 
a  deep  tunnel  (500m  deep),  and  3)  FE  models  of  the  SRI  SWAT  precision 
tests,  and  4)  a  shallow  tunnel  targeting/vulnerability  problem.  In  these 
analyses,  the  advanced  probabilistic  methods  were  shown  to  be  accurate 
and  orders  of  magnitude  more  efficient  than  Monte  Carlo  simulation. 

Many  of  the  algorithms  enhanced  and/or  implemented  on  the  DNA  project  were  developed  to 
some  degree  by  a  9-year  NASA  project  with  SwRI  entitled  "Probabilistic  Structural  Analysis 
Methods  for  Select  Space  Shuttle  Components"  (PSAM).  PSAM  was  aimed  at  developing 
probabilistic  strucmral  analysis  methods  for  space  shuttle  main  engine  applications.  This 
significant  level  of  technology  transfer  from  this  NASA  project  to  the  DNA  project  greatly 
increased  the  scope  and  magnitude  of  accomplishments  than  what  would  have  otherwise  been 
possible. 
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SECTION  1 


INTRODUCTION 


1.1  BACKGROUND. 

Predicting  the  behavior  and  failure  of  deep  underground  tunnels  in  severe  loading  environments 
is  now  becoming  possible  with  new  computational  techniques  and  more  powerful  computers. 
The  problems  are  highly  transient,  involve  nonlinear  material  behavior  and  include  large 
uncertainties  in  material  properties,  loading,  and  geometry.  Sophisticated  numerical  models 
based  on  finite  and  distinct  element  methods  are  being  developed  and  applied  to  these  problems. 
However,  computational  costs  are  high  as  even  a  single  analysis  may  require  hours  of  computer 
time.  This  computational  cost  coupled  with  the  need  to  account  for  uncertainties  in  a  direct 
manner  has  led  to  the  development  and  application  of  efficient  probabilistic  analysis  methods. 

Parameters  used  to  describe  geomechanics  problems  are  known  to  possess  large  uncertainties. 
These  large  uncertainties  in  material  properties,  loading,  and  geometry  stem  from  the  fact  that 
the  composition  of  most  rock  is  highly  variable  and  includes  joints  and  faults  not  accounted  for 
in  material  testing.  An  assessment  of  the  variability  or  randomness  of  the  structural  behavior 
and/or  failure  requires  careful  study  for  a  thorough  understanding  of  the  problem.  The 
traditional  method  of  accounting  for  uncertainties  is  a  factor-of-safety  approach.  However,  this 
approach  does  not  quantify  the  failure  probability  of  a  given  design,  which  can  lead  to  either 
overconservative  or  unreliable  designs,  nor  does  it  identify  the  variable  contributing  the  most  to 
the  failure  probability. 

The  traditional  method  of  probabilistic  analysis  is  Monte  Carlo  simulation.  This  approach 
generally  requires  hundreds  of  thousands  of  simulations,  and  is  impractical  when  each 
simulation  involves  extensive  numerical  calculations.  The  pursuit  of  alternative  probabilistic 
analysis  methods  that  are  more  efficient  than  Monte  Carlo  simulation  has  led  to  the  development 
of  many  new  approximate  probabilistic  methods  over  the  past  two  decades.  Many  of  these 
methods  have  been  shown  to  be  highly-accurate  and  applicable  for  use  with  complex  numerical 
deterministic  models. 

1.2  OBJECTIVE. 

The  objective  of  this  research  was  to  develop  a  methodology  that  predicts  the  uncertainty  in 
structural  deformation  and  damage  in  underground  tunnels  given  the  uncertainty  in  geologic  and 
reinforcing-structure  parameters.  The  method  was  to  be  more  efficient  than  traditional  Monte 
Carlo.  The  development  of  fast-running  engineering  models  for  tunnel  response  was  also 
required. 
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The  approach  taken  was  to: 

1 .  Develop  a  three-dimensional  numerical  model  for  the  response  of  deep-buried  tunnels. 

2.  Verify  and  validate  the  deterministic  numerical  model. 

3.  Develop  and  integrate  efficient  probabilistic  methods  with  the  validated  deterministic 
model. 

4.  Prepare  probabilistic  inputs. 

5.  Calculate  the  uncertainty  in  structural  deformation  and  identify  key  problem  variables  for 
several  selected  problems. 

6.  Demonstrate  probabilistic  analysis  methods  using  a  fast-running  tunnel  targeting  model. 
1.3  SCOPE. 

To  satisfy  the  objective,  several  new  and  efficient  probabilistic  methods  were  developed  and 
integrated  with  a  general  purpose  numerical  model.  A  common  element  to  the  methods  is  the 
use  of  fast  probability  integration  (FPI)  algorithms,  which  have  been  shown  to  be  many  times 
more  efficient  than  Monte  Carlo  simulation  (Wu  and  Wirsching,  1987).  The  advanced  mean- 
value  (AMV)  procedure,  which  is  based  on  FPI,  was  found  to  be  extremely  efficient  in 
computing  the  probabilistic  response  of  complex  structures  with  relatively  few  response 
calculations.  As  a  means  of  verification,  several  advanced  simulation  methods  were  also 
employed  such  as  Latin  Hypercube  Simulation  (LHS)  and  Adaptive  Importance  Sampling  (AIS) 

This  report  is  organized  as  follows:  Section  2  summarizes  the  development  of  two  deterministic 
models;  an  analytical  tunnel  closure  model  and  a  more  general  purpose  finite  element  (FE) 
model.  Section  3  reports  on  the  verification  and  validation  of  the  FE  model.  In  Section  4,  the 
probabilistic  methods  employed  are  presented.  The  probabilistic  FE  model  is  demonstrated  on 
several  problems  in  Section  5,  including  a  deep  tunnel  analysis,  a  dynamic  laboratory 
experiment,  and  a  tunnel  vulnerability  problem.  References  to  the  literature  are  used  where 
possible  to  keep  the  report  concise. 
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SECTION  2 


DETERMINISTIC  TUNNEL  MODELS 


2.1  INTRODUCTION. 

In  this  Section,  two  deterministic  models  are  presented.  The  first  model,  based  on  the  work 
of  Hendron-Aiyer  (1972),  has  been  used  for  tunnel  analysis  and  includes  both  elastic  and 
plastic  closure.  Although  the  model  makes  several  simplifying  assumptions,  it  has  been  used 
extensively  in  the  DNA  community  for  making  rapid  assessments  and  has  performed 
remarkably  well.  The  second  model  is  a  general-purpose  finite  element  solid  dynamics 
program  (PRONTO)  employing  several  constitutive  models  for  modeling  the  elastic  and 
plastic  response  of  intact  rock  and  joints.  Since  PRONTO  is  developed  and  maintained  by 
Sandia  National  Laboratories  (SNL)  and  is  available  for  government  and  commercial  use,  the 
enhancements  made  to  PRONTO  during  the  course  of  this  program  are  readily  available. 

2.2  ANALYTICAL  ELASTOPLASTIC  TUNNEL  CLOSURE  MODEL. 

I 

The  problem  solved  is  that  of  a  long,  circular  tunnel  of  radius  a  in  an  elastic-perfectly  plastic 
Mohr-Coulomb  material  subjected  to  static  axisymmetric  pressure  on  the  boundary  of  the 
tunnel  and  also  at  infinity.  Plane  strain  is  assumed  along  the  axis  of  the  tunnel,  and  the 
plastic  flow  of  the  medium  may  be  nonassociative.  The  solution  assumes  that  the  internal 
and  far-field  pressure  are  increased  together  until  the  desired  internal  pressure  is  imposed. 

The  far-field  pressure  is  then  further  increased,  resulting  in  yielding  of  the  medium  around  the 
tunnel  and  increased  tunnel  closure. 

The  problem  geometry  is  given  in  Figure  2-1.  The  variables  in  the  model  are  Young's 
modulus,  Poisson's  ratio,  unconfined  compressive  strength,  angle  of  internal  fnction,  a  dilation 
parameter,  internal  pressure,  and  external  pressure. 

The  formulae  for  tunnel  closure,  AD/D,  are  given  in  Thacker  and  Senseny  (1992)  from  the 
solution  for  the  stress  and  displacement  fields  obtained  by  Wintergerst,  et  ah,  (1991).  Their 
work  extended  the  solutions  obtained  by  Florence  and  Schwer  (1978)  to  allow  for  non- 
associated  plastic  flow.  The  Florence  and  Schwer  solution  removed  the  constraint  of  earlier 
solutions  (Newmark,  1969;  Hendron  and  Aiyer,  1972)  that  the  least  compressive  principal 
stress  be  the  radial  stress.  A  complete  description  of  the  model  and  the  analytical  tunnel 
closure  equations  can  be  found  in  Thacker  and  Senseny  (1992). 
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(a)  -  Outer  Elastic  and  Three  Plastic  Zones  (b)  -  Outer  Elastic  and  Two  Plastic  Zones 

Figure  2-1.  Analytical  tunnel  problem  geometry:  a)  case  1  -  outer  elastic  and  three  plastic  zones,  b) 
case  2  -  outer  elastic  and  two  plastic  zones. 


2.3  EXPLICIT  DYNAMIC  FINITE  ELEMENT  MODEL. 

PR0NT02D  is  a  two-dimensional  transient  solid  dynamics  code  for  analyzing  large 
deformations  of  highly  nonlinear  materials  subjected  to  extremely  high  strain  rates  (Taylor  and 
Flanagan,  1987).  It  is  the  latest  in  a  series  of  transient  dynamic  finite  element  codes  that  have 
been  developed  at  Sandia  National  Laboratories,  beginning  with  HONDO  (Key,  et  ai,  1978). 

As  such,  PR0NT02D  contains  a  number  of  state-of-the-art  numerical  algorithms,  including  an 
adaptive  timestep  control  algorithm,  a  robust  hourglass  control  algorithm,  a  very  accurate 
incremental  rotation  algorithm,  and  a  robust  surface  contact  algorithm.  Four-noded,  uniform- 
strain,  quadrilateral  elements  with  single-point  integration  are  used  in  the  finite  element 
formulation.  Beyond  its  general  capabilities,  PR0NT02D  was  chosen  for  the  calculations 
because  new  constitutive  models  are  readily  added  and  because  a  three-dimensional  derivative  of 
the  program  is  available  (PR0NT03D). 

The  two  features  that  make  the  deep  tunnel  calculations  relatively  unique  among 
large-deformation  solid  dynamics  problems  are:  1 )  a  jointed  rock  mass  subject  to  loading 
conditions  that  could  result  in  large  relative  motions  (sliding)  between  adjacent  rock  blocks,  and 
2)  loading  conditions  that  could  result  in  substantial  tensile  failure,  plastic  deformation,  and 
dilation  within  the  rock  blocks.  The  joints  also  affect  the  elastic  behavior  of  the  rock  mass  by 
reducing  its  elastic  moduli  from  the  intact  rock  values.  Further,  the  elastic  moduli  of  jointed 
rock  generally  are  nonlinear  functions  of  the  joint  apertures. 
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The  surface  contact  algorithm  in  PRONT02D  is  directly  applicable  to  modeling  the  slippage 
between  adjacent  rock  blocks.  It  is  designed  to  simulate  smooth,  cohesionless  surfaces  with 
shear  strengths  defined  in  terms  of  static  and  dynamic  coefficients  of  friction.  In  the  tunnel 
problems,  the  static  coefficient  of  friction  is  the  tangent  of  the  joint  friction  angle.  The  friction 
coefficient  does  not  decay  during  sliding,  so  the  dynamic  coefficient  of  friction  is  equal  to  the 
static  coefficient. 


The  surface  contact  algorithm  does  not  provide  a  means  to  model  the  elastic  behavior  of  the 
joints.  Consequently,  the  approach  used  is  to  superpose  the  nonlinear  stress-displacement 
response  of  the  joints  and  the  linear  stress-strain  response  of  the  intact  rock.  The  result  is  a 
composite  material  whose  elastic  behavior  is  equivalent  to  a  jointed  rock  mass.  This  approach 
has  been  developed  and  used  successfully  by  other  modelers  (Labreche  and  Petney,  1987)  and  it 
is  referred  to  as  the  Complaint  Joint  Model  (CJM).  Note  that  the  elastic  behavior  of  jointed  rock 
masses  and  as  simulated  by  the  CJM  is  inherently  anisotropic  because  the  stress-strain  response 
normal  to  a  joint  set  is  substantially  different  from  the  response  tangential  to  the  joint  set. 


2.3.1  Explicit  Time  Integration. 


The  equations  of  motion  are  integrated  using  a  modified  central  difference  scheme  in 
PR0NT02D.  The  velocities  are  integrated  with  a  forward  difference  and  the  displacements  are 
integrated  with  a  backward  difference.  This  scheme  is  expressed  as 


= 


^t^Li  ^  DtT 

Jt  ^  Jt 
m 


+  Ar 


(2.1) 


+  =  “r  + 


A/U 


i  +  Ar 


where  is  the  external  nodal  force,  ^  is  the  internal  nodal  force,  m  is  the  nodal  point 
lumped  mass,  and  Af  is  the  time  increment.  This  central  difference  operator  is  conditionally 
stable  and  the  Courant  stability  limit  is  given  by  the  highest  eigenvalue  of  the  system 
(Bathe  and  Wilson,  1976) 


Af  i 


(2.2) 


Flanagan  and  Belytschko  (1984)  provided  eigenvalue  estimates  for  the  uniform  strain 
quadrilateral  used  in  PRONT02D. 

Numerical  damping  is  introduced  in  the  solution  by  adding  artificial  viscosity.  This  prevents 
high  velocity  gradients  from  collapsing  an  element  before  is  has  a  chance  to  respond  and  to  quiet 
truncation  frequency  ringing.  The  technique  used  in  PR0NT02D  is  to  add  viscosity  to  the 
"bulk"  response.  This  generates  a  viscous  pressure  in  terms  of  the  volume  strain  rate. 
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2.3.2  Four  Node  Uniform  Strain  Element. 


PR0NT02D  uses  a  four-noded  two-dimensional  uniform  strain  element  in  the  finite  element 
formulation.  A  one  point  integration  of  the  element  under-integrates  the  element  but  provides  a 
large  computational  advantage  over  a  two-by-two  integration  rule.  However,  this  results  in  a 
rank  deficiency  for  the  element  that  may  cause  spurious  zero  energy  (hourglass)  modes. 

The  uniform  strain  element  considers  only  a  fully  linear  field.  Any  remaining  nodal  velocity 
field  is  the  hourglass  field.  Possibly  severe  unrestricted  mesh  distortion  can  occur  if  these  modes 
are  excited.  The  method  used  in  PR0NT02D  isolates  the  hourglass  modes  so  they  may  be 
treated  independently  of  the  rigid  body  and  uniform  strain  modes  (Flanagan  and  Belytschko, 
1981).  PR0NT03D,  which  also  uses  reduced  integration  elements,  uses  a  similar  approach  to 
control  hourglassing. 

2.3.3  Material  Behavior. 

2.3.3. 1  Dnicker-Prager.  Several  classical  yield  surfaces  defined  in  terms  of  the  first  stress 

invariant  and  the  second  deviatoric  stress  invariant  have  been  used  to  model  plasticity  in  rock. 
They  include  the  linear  Mohr-Coulomb  and  the  Drucker-Prager  yield  criteria.  The  Drucker- 
Prager  yield  function  is 

f(o^)  =  /jl-aJ.-k  (2.3) 


where  o  are  the  components  of  the  stress  tensor,  a  and  k  are  material  constants,  is  the  first 
invariant  of  the  stress  tensor,  and  second  invariant  of  the  deviatoric  stress  tensor.  The 

yield  surface  is  the  locus  of  stress  states  at  which  the  value  of  the  yield  function  is  zero  (/=  0). 
Consequently,  the  Drucker-Prager  yield  surface  is  a  cone  in  principal  stress  space  with  the  axis 
of  the  cone  along  the  hydrostat. 


The  Drucker-Prager  and  the  Mohr-Coulomb  yield  criteria  are  equivalent  only  at  certain  stress 
states,  depending  on  how  the  Drucker-Prager  material  constants  are  evaluated.  To  match  the 
Mohr-Coulomb  criterion  in  triaxial  compression,  the  material  constants  must  be  evaluated  as 
follows: 


_  2sm<() 

V^(3-sm4)) 


it  = 


2v^cos<{)  ^ 
3-sin4)  j  ® 


(2.4) 


where  (j>  andc^  are  the  Mohr-Coulomb  friction  angle  and  cohesion,  respectively.  The  resultant 
Drucker-Prager  yield  surface  circumscribes  the  Mohr-Coulomb  yield  surface. 

The  stress  state  is  elastic  when  the  yield  function  is  negative  (/<  0).  When  the  stress  state  is  on 
the  yield  surface  (/=  0)  and  the  loading  condition  is  such  that  it  would  result  in  /iO,  the 
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resulting  deformation  is  plastic.  The  components  of  the  resultant  plastic  strain  tensor  are  defined 
in  terms  of  the  flow  rule,  which  is  classically  stated  as  proportional  to  the  gradient  of  a  plastic 
potential  function  g(o^).  The  proportionality  constant  is  determined  by  the  consistency 
condition  which  requires  that  the  stress  state  remain  on  the  yield  surface  during  plastic 
deformation.  When  the  plastic  potential  function  and  the  yield  function  coincide  for  all  stress 
states  (g  =f)  the  flow  rule  is  said  to  be  associated  with  the  yield  function. 


An  incremental  method  using  tangent  stiffness  (Chen  and  Han,  1988)  is  used  to  implement 
Drucker-Prager  plasticity  in  the  constitutive  model  that  was  added  to  PR0NT02D  for  the 
benchmark  problems.  TTie  incremental  stress-strain  relationship  can  be  expressed  in  the 
following  general  form: 


da„  = 


dg  df 


ao,. 

df  _  a* 


da. 


rstu 


da 


tu 


dsu 


(2.5) 


where  da^  and  are  the  components  of  the  incremental  stress  and  strain  tensors,  respectively, 
are  the  components  of  the  elastic  coefficient  tensor,  and  repeated  subscripts  indicate 
summation.  The  coefficient  tensor  in  the  brackets  represents  the  elastic-plastic  tensor  of  tangent 
moduli  for  an  elastic-perfectly  plastic  material.  The  quantities  in  the  brackets,  including  for 
a  nonlinear  elastic  material  such  as  in  the  CJM,  are  evaluated  at  the  current  stress  state. 


For  the  Drucker-Prager  yield  surface,  the  partial  derivatives  of  the  yield  function  with  respect  to 
the  stress  components  are 


df 


(2.6) 


where  are  the  components  of  the  deviatoric  stress  tensor,  (Sj  =  a^-  6„/j/3)and  6^  is  the 
Kronecker  delta  function.  In  the  constitutive  model  added  to  PR0NT02D  the  following 
extended  form  of  the  plastic  potential  function  associated  with  the  Drucker-Prager  yield  function 
is  used: 


where 


gio^)  = 


2sinT|f 

/J(3  -smi|f) 


(2.7) 


(2.8) 
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and  i|r  is  the  dilation  angle.  The  plastic  potential  function  results  in  an  associated  flow  rule 
when  the  dilation  angle  equals  the  friction  angle  (so  that  p  =  a  and  g  =/).  The  components  of 
the  flow  rule  (the  partial  derivatives  of  g  with  respect  to  a^)  are 


(2.9) 


2.3.3. 2  Cap  Model.  The  Sandler-Rubin  (SR)  cap  model  proposed  in  Sandler  and  Rubin  (1979) 
is  an  isotropic  rate-independent  elastic-plastic  material  model  based  on  the  classical  incremental 
theory  of  plasticity.  It  has  been  used  to  model  both  high  and  low  pressure  mechanical  behavior 
of  many  geologic  materials  such  as  sands,  clays  and  rocks  because  of  its  ability  to  simulate  the 
important  characteristics  of  these  materials:  compaction,  dilatancy,  shear  and  hysteresis.  The 
yield  surface  consists  of  a  fixed  failure  envelope  given  by 


=  /^'  =  A-Cexp[B'J^\  for  /j  >  L(k)  (2.10) 

t 

and  a  hardening  yield  function  (or  cap)  given  by 

Fj,(7i.k)  =  /^'  =  ^^[X(K)-L(K)f-[/,-L(K)f  for  L(K)i/iiX(K)  (2.11) 


where  A,  B,  and  C  are  the  parameters  defining  the  shear  failure  surface,  /j  is  the  first  invariant 
of  the  stress  tensor,  is  the  second  invariant  of  the  deviatoric  stress  tensor,  k  is  a  hardening 
variable  related  to  the  total  plastic  volumetric  strain,  X(k)  and  L(k)  are  the  values  of  at  the 
intersection  of  the  cap  with  the  /,  axis  and  center  of  the  cap  respectively,  and  R  is  the  ratio  of 
the  major  to  minor  axes  of  the  elliptical  cap.  In  (2.1 1), 


if  K  <  0 
1/  K  i  0 


(2.12) 


and 

X(k)  =  LiK)-R{A-Cexp[B'L(K)]}  =  L(k)-R-Fj,(L(k)).  (2.13) 

The  hardening  variable  k  is  related  to  the  plastic  volume  strain  (caused  only  by  cap  action) 
through  the  relationship 

^  =  W{cxp  [D  (X(K)  -X;,)]  - 1}  (2.14) 
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Where  XJ,  is  the  initial  location  of  the  cap,  and  to  keep  the  cap  from  retracting  (for  rocks) 

(2.15) 


The  original  SR  cap  was  generalized  to  accommodate  the  anisotropic  elastic  CJM  model. 


Figure  2-2.  Sandler-Rubin  cap  model. 


The  failure  surface  of  the  SR  cap  model  is  also  modified  to  account  for  different  strengths  in  triaxial 
extension  (TXE)  and  triaxial  compression  (TXC).  Specifically,  the  failure  surface  defined  by 

f  =  ^2-  ai-r  (2.16) 


is  modified  by  31  resulting  in  a  smoothly-varying  yield  surface  in  the  it -plane,  as  shown  in 
Figure  2-3.  31  is  given  by 

31  =  -[1  -»-r-l-(l  -  X)sin30]  (2.17) 

2 


where  0  is  the  Lode  angle  which  is  related  to  the  third-invariant  of  the  deviatoric  stress  tensor, 
and  X=TXE/TXC  limited  to  0.778 1.0  to  ensure  a  convex  yield  surface.  Also  shown  in 
Figure  2-3  are  the  familiar  Mises-Schleicher  (MS)  and  Mohr-Coulomb  (MC)  failure  surfaces. 
Note  that  the  MS  surface  is  obtained  by  fitting  the  model  to  TXC  conditions  and  does  not 
account  for  any  reduced  strength  in  TXE.  MC,  on  the  other  hand,  is  defined  by  straight  lines 
between  TXE  and  TXC.  Because  some  of  the  calculations  performed  by  other  DNA  contractors 
(presented  in  Section  3)  were  performed  using  the  Willam-Wamke  (1975)  function,  it  is  also 
shown  in  Figure  2-3  to  illustrate  the  differences  between  it  and  the  failure  surfaces  just 
described. 
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It  is  well  known  that  materials  exhibit  higher  strengths  when  loaded  at  higher  strain  rates. 

Because  shock-loaded  tunnels  involve  strain-rates  as  high  as  10  sec'  ,  rate  effects  must  be 

accounted  for.  Unfortunately,  only  a  limited  amount  of  data  is  available  at  these  high 

strain-rates.  Green  and  Perkins  (1972)  and  others  have  measured  a  sharp  increase  (exponential) 

2  “1 

in  strength  for  strain  rates  above  about  10  sec  . 

In  the  high  compressive  states  associated  with  shock  loadings,  rock  becomes  an  almost  plastic 
material.  Given  the  success  of  rate-independent  cap  models  for  representing  this  behavior,  a 
number  of  researchers  have  used  two  approaches  in  developing  viscoplastic  versions  of  the  cap 
model  to  represent  rock  under  high  loading  rates.  In  the  approach  of  Perzyna  (1966),  the 
viscoplastic  strain  rate  is  a  function  of  the  level  of  overstress  and  a  fluidity  parameter  A. , 

=  ^  ((!>(/))  ^  (2.18) 

do„ 


In  (2.18),  4>(/)  can  be  chosen  to  match  experimental  data.  In  the  Duvaut  and  Lions  (1972) 
approach,  the  viscoplastic  strain  rate  is  a  linear  function  of  the  overstress,  and  a  material  constant 
termed  the  relaxation  time  x  that  determines  the  rate  at  which  the  stress  relaxes  back  to  the 
inviscid  solution.  The  Duvaut-Lions  model  is  given  by 

4*  =  (2.19) 

X 
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where  is  the  elasticity  tensor,  is  the  stress  tensor  and  is  the  inviscid  stress  tensor.  The 
Duvaut-Lions  model  is  attractive  because  of  its  simplicity  and  ease  of  implementation. 

2.3.3.3  T.imited  Tension  Algorithm.  A  limited-tension  algorithm  developed  by  Callahan  (1989) 
is  used  in  the  constitutive  model  to  simulate  tensile  failure.  In  this  algorithm  the  tensile  strength 
at  each  integration  point  is  set  initially  to  the  intact  tensile  strength  At  each  integration 
point  at  each  time  step,  the  principal  stresses  o^)  and  their  unit  direction  vectors 

(Aj,  ftjj,  ijjf)  are  calculated  based  on  the  trial  stress  state  oj  at  that  location  and  time.  Each 
principal  stress  component  that  exceeds  the  tensile  strength  at  the  integration  point  is  set  to 
the  residual  tensile  strength  T^.  In  addition,  the  tensile  strength  at  the  integration  point  is  set  to 
if  r,  is  exceeded.  (Recall  that  tension  is  negative  according  to  the  sign  convention,  so 
exceeding  implies  that  the  stress  is  less  than  T,  and  the  stress  is  set  to  -T^.)  The  resultant 
stress  state  is  calculated  by  transforming  the  principal  stress  components  back  to  the  global 
coordinate  system  while  maintaining  the  original  principal  stress  directions.  This  transformation 
is  accomplished  according  to  the  following  equation: 


(2.20) 


where  is  the  i*  component  of  unit  direction  vector  Aj. 

2,3.3.4  Implicit  Rock  Joint  Modelling.  The  compliant  joint  model  (CJM)  is  a  three-dimensional 
elastic  model  used  to  represent  the  mechanical  response  of  a  rock  matrix  with  joints  (Labreche 
and  Pemey,  1987).  The  model  implemented  in  PR0NT02D  is  a  two-dimensional  version  of  the 
CJM.  The  CJM  assumes  that  the  deformation  of  a  jointed  rock  mass  is  caused  by  the 
combination  of  the  elastic  response  of  the  rock  matrix  and  the  elastic  response  of  the  joints.  The 
unfractured  rock  matrix  is  modeled  as  an  isotropic  linear  elastic  solid  and  the  joints  as  a 
nonlinear  elastic  layer.  The  model  consists  of  up  to  four  joint  sets  with  various  spacing  and 
orientation.  The  response  of  the  CJM  is  effectively  anisotropic  due  to  the  presence  of  the  joint 
sets. 

The  displacement  behavior  normal  to  the  joint  plane  varies  nonlinearly  with  stress  and  the 
response  is  governed  by  the  physical  characteristics  such  as  joint  aperture,  roughness  and  contact 
between  joint  surfaces,  strength  of  the  wall  rock,  and  presence  or  absence  of  filling.  The 
response  is  also  governed  by  the  recent  movement  history  of  the  rock  joint.  The  aperture  change 
of  a  well-mated  joint  without  filling  can  be  represented  by  a  hyperbolic  relationship  between  the 
stress  normal  to  joint  plane  and  the  joint  closure  (Goodman,  1976;  Bandis,  et  al.,  1983).  The 
form  used  in  the  CJM  is  given  by  the  equation 


o 


a 


/ 

a 

\ 


6  +  V 


■/ 


(2.21) 


where, 

o  -  normal  stress  acting  on  the  joint  (compression  positive) 
a  -  half-closure  stress  (stress  to  reduce  the  aperture  to  6  /2 


11 


V  -  joint  normal  displacement  (closure  positive,  zero  when  o,  0) 
6  -  maximum  joint  closure 

m  -  exponent 


The  stiffness  normal  to  the  joint,  k^,  is  defined  by 

k  =  — -  =  am  - 


(2.22) 


which  is  the  slope  of  the  nonlinear  joint  closure.  This  hyperbolic  function  is  used  in  the 
numerical  formulation  for  the  normal  stress-displacement  response  when  the  normal  stress  is 
compressive.  A  constant  value  of  10  ^  times  the  dilatational  modulus,  X  ~2ii,  of  the  matrix  is 
used  for  joint  normal  stiffness  when  the  normal  stress  is  tensile.  This  stiffness  defines  a  linear, 
normal  stress-normal  displacement  relationship  for  the  joint  in  tension.  This  small  but  nonzero 
stiffness  avoids  numerical  singularities  in  the  computations.  The  transmission  of  small  tensile 
stresses  by  elements  having  moduli  four  orders  of  magnitude  lower  than  the  matrix  are  accepted 
as  insignificant  relative  to  typical  compressive  stresses. 

The  compressive  normal  stress-displacement  relation  is  nonlinear  and  elastic.  Thus,  there  is  no 
path  dependence  and  all  joint  closure  is  fully  recoverable  in  this  model. 

The  shear  response  of  joints  generally  consists  of  a  nearly  linear  rise  to  peak  shear  strength, 
followed  by  softening  behavior  with  continuing  shear.  However,  it  is  acknowledged  that  joint 
size  dependence  may  alter  this.  The  peak  shear  strength  response  of  joints  is  supported  by 
studies  such  as  Goodman  (1976)  and  Bandis  (1983).  The  shear  stress  of  joints  is  known  to  be 
inelastic  when  the  stress  state  approaches  the  peak  shear  strength.  The  CJM  only  models  the 
elastic  joint  shear  response.  The  shear  stiffness,  which  is  independent  of  the  normal  stiffness 
across  the  joint,  relates  the  shear  displacement  to  shear  stress. 

X  =  k,V^  (2.23) 


where. 


X  -  shear  strength 

k^  -  joint  shear  stiffness 

Vj  -  joint  displacement  tangent  to  plane  of  joint 

The  shear  stiffness  is  defined  by 


(2.24) 
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The  CJM  introduces  anisotropy  by  reducing  the  composite  modulus  of  the  modeled  rock  mass  in 
the  direction  normal  to  the  joint  in  an  otherwise  linear  elastic  material.  Computation  of  effective 
moduli  is  used  to  implement  the  softening  effects  of  the  joints.  Stress  equilibrium  is  the 
underlying  assumption  in  the  combined  response  of  the  rock  matrix  and  the  joint  sets.  The 
normal  stress  across  the  joint  is  equal  to  the  stress  in  the  rock  matrix  normal  to  the  joint  plane. 
The  shear  stress  on  the  joint  is  equal  to  the  shear  stress  in  the  rock  matrix  parallel  to  the  joint 
plane. 

Correspondingly,  the  strains  in  the  rock  matrix  are  combined  with  the  strains  in  the  joints 
additively.  That  is,  the  effective  compliance  for  the  rock  mass  is  the  sum  of  the  compliances  of 
the  unfractured  rock  and  the  joint  sets.  The  effective  Young's  modulus,  t,  normal  to  a  joint  and 
the  effective  shear  modulus,  6,  parallel  to  a  joint,  are 

1  =  1  +  J. 

^  £  V 

(2.25) 

A  =  _L  +  J_ 

d  G  kjs 

where  E  and  G  are  the  elastic  moduli  of  the  rock  matrix,  and  k^  are  the  normal  and  shear 
stiffness  of  the  joint  set,  and  s  is  the  spacing  of  the  joints  in  the  joint  set. 

2.3.3. 5  Explicit  Rock  Joint  Modelling.  PRONTO  treats  contact  as  a  kinematic  constraint.  That 
is,  the  algorithm  modifies  the  accelerations  of  the  nodes  along  the  contact  region  such  that  the 
kinematic  constraints  are  satisfied.  The  algorithm  uses  a  partitioned  approach  to  enforce 
compliance  between  two  contact  surfaces  by  allowing  each  surface  to  act  as  a  master  surface  for 
a  fraction  of  each  time  step  and  a  slave  for  the  remainder. 

There  are  four  steps  involved  in  the  contact  algorithm.  First,  the  contact  surface  geometry  is 
recalculated  at  each  time  step  and  the  predicted  configuration  is  computed  by  integrating  the 
motion  without  regard  to  the  kinematic  constraints  required  by  the  contact  surfaces.  The 
following  quantities  are  calculated  for  each  node: 

d  =  //m 

«  =  v+AM  (2.26) 

t  =  x  +  AfO 


Where  /  is  the  residual  force  vector  (sum  of  external  forces  minus  the  sum  of  internal  forces),  m 
is  the  nodal  mass,  v ,  is  the  current  velocity,  and  Af  is  the  time  increment.  The  predicted 
kinematic  quantities  are  denoted  by  the  hat. 
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The  second  step  is  surface  tracking,  or  the  process  of  matching  nodes  along  one  surface  with  the 
mating  surface.  The  algorithm  used  is  to  locate  the  spatially  nearest  master  node  to  the  possible 
point  of  contact.  This  procedure  can  be  the  most  time  consuming  portion  of  the  analysis  for  this 
type  of  problem.  To  streamline  the  tracking  algorithm,  the  nearest  master  node  to  a  given  slave 
node  at  one  time  step  is  assumed  to  be  in  the  vicinity  of  the  nearest  master  node  at  the  next  time 
step.  Therefore,  at  each  time  step,  the  nearest  master  node  is  the  starting  point  for  the  search 
along  the  surface. 

The  next  step  is  to  determine  contact  or  penetration.  The  slave  node  is  oriented  with  respect  to 
the  master  segments  connected  to  the  tracked  master  node.  The  depth  and  position  coordinates 
are  calculated  for  both  master  segments  connected  to  the  nearest  master  node.  From  these 
quantities,  if  the  depth  is  positive,  the  slave  node  is  penetrating  the  segment  and  if  the  position  is 
positive  then  the  slave  node  is  along  the  segment.  When  the  master  surface  forms  an  outside 
comer,  there  may  be  penetration  of  both  segments.  In  this  case,  the  algorithm  determines  with 
which  master  segment  the  slave  surface  is  more  strongly  in  contact.  One  limitation  of  the 
algorithm  is  that  it  can  not  detect  a  contact  surface  contacting  itself. 

The  final  step  in  the  contact  algorithm  is  to  calculate  the  penetration  forces  imposed  on  the 
master  surface  by  the  slave  surface  to  restore  kinematic  compliance.  These  forces  are  calculated 
as  a  fraction  of  the  forces  that  would  be  imposed  by  the  slave  nodes  if  the  master  surface  was 
rigid.  This  fraction  is  based  on  the  fraction  of  each  time  step  for  which  the  surfaces  act  as 
master  and  slave.  The  roles  are  reversed  for  the  remaining  fraction  of  the  time  step.  The 
accelerations  are  calculated  to  predict  the  response  of  the  master  surface  to  these  penetration 
forces  such  that  the  response  of  each  contacting  slave  node  is  constrained  by  its  master  nodes. 
The  principle  of  virtual  work  is  employed  to  define  the  accelerations  of  the  master  nodes  in 
response  to  the  penetration  forces.  When  friction  is  present,  the  relative  tangential  motion  of  the 
contacting  slave  nodes  is  resisted.  The  magnitude  of  the  tangential  force  exerted  on  the  master 
surface  on  a  slave  node  cannot  exceed  the  friction  force.  Thus  friction  adds  a  tangential 
acceleration  to  the  nodes  in  contact. 

2.3.4  Transmitting  Boundary  Conditions. 

A  transmitting  boundary  is  used  to  simulate  a  semi-infinite  domain  outside  the  boundary,  where 
the  wave  speeds  of  the  material  on  both  sides  of  the  boundary  are  the  same.  The  region  exterior 
to  the  boundary  is  replaced  with  an  energy-absorbing  boundary  condition  that  behaves  as  if  the 
energy  is  transmitted  across  the  boundary.  Thus,  no  energy  is  reflected  back  into  the  interior 
region. 

The  nonreflecting  boundary  is  implemented  in  PRONTO  according  to  a  technique  proposed  by 
Lysmer  (Lysmer  and  Kuhlemeyer,  1979).  The  basic  idea  is  to  apply  tractions  at  the  boundary 
that  will  exactly  cancel  the  stresses  that  would  be  reflected  from  a  free  surface.  Hence,  the 
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numerical  technique  involves  calculating  and  applying  the  following  tractions  to  the  surface  of 
the  nonreflecting  boundary 


^  =  py.K 

(2.27) 


where  and  are  the  normal  and  shear  tractions;  p ,  Vp ,  and  are  the  current  density, 
compressional  wave  speed,  and  shear  wave  speed  of  the  material  along  the  boundary;  and  ti^  and 
il,  are  the  normal  and  tangential  components  of  the  current  velocity  at  the  boundary.  At  each 
time  step,  PRONTO  updates  the  tractions  at  the  boundary  using  the  current  density  and  effective 
dilatational  and  shear  moduli  in  each  element  along  the  boundary. 

In  calculating  the  current  wave  speeds,  PRONTO  assumes  that  the  wave  speeds  are  isotropic  and 
independent  of  the  orientation  of  the  boundary.  Obviously,  the  joint  sets  in  a  rock  mass  yield  an 
anisotropic  medium  in  which  the  wave  speeds  depend  on  the  incident  direction.  However,  the 
anisotropy  inherent  in  the  CJM  is  not  accounted  for  in  the  wave  speed  calculation  in  PRONTO. 
Consequently,  the  wave  speeds  used  to  calculate  the  tractions  along  the  nonreflecting  boundary 
are  somewhat  in  error,  the  magnitude  of  which  depends  on  the  orientation  of  the  boundary  with 
respect  to  the  principal  directions  of  the  anisotropy.  This  error  results  in  a  partial  reflection  of 
the  incident  stress. 


15 


SECTION  3 


VERIFICATION  AND  VALIDATION  OF  DETERMINISTIC  TUNNEL  MODELS 


3.1  INTRODUCTION. 

The  quality  of  probabilistic  predictions  is  highly  dependent  on  the  quality  of  the  underlying 
deterministic  model.  For  complex  problems  that  have  no  analytical  solution,  the  quality  of  the 
deterministic  model  can  be  measured  qualitatively  in  terms  of  confidence  gained  by  applying  the 
model  and  comparing  predicted  results  with  test  measurements  for  a  variety  of  different 
configurations  and  loadings. 

The  process  of  obtaining  confidence  in  the  deterministic  model  can  be  termed  verification  and 
validation.  Verification  involves  comparing  predictions  to  either  known  solutions,  or  comparing 
predictions  that  have  no  known  solution  to  independently  made  predictions  using  other  methods. 
While  the  former  is  preferable,  comparing  to  independently  obtained  solutions  can  be  an 
effective  way  to  augment  the  suite  of  verification  problems.  \Vhen  a  numerical  model  is  said  to 
be  verified,  it  is  assumed  to  be  working  as  intended  and  be  free  firom  errors,  bugs,  and  numerical 
inaccuracies. 

The  process  of  validating  a  numerical  model  involves  comparing  predictions  made  using  the 
numerical  model  with  experimental  results.  The  process  is  much  more  difficult  than 
verification,  and  in  reality,  can  never  be  achieved  completely.  However,  a  model  is  usually 
accepted  as  being  validated  if  it  has  been  demonstrated  to  be  able  to  predict  all  significant 
physical  processes  observed  and  measured  in  a  similar  test. 

To  establish  confidence  in  the  numerical  models,  SwRI  participated  in  two  large  verification  and 
validation  programs,  which  are  sununarized  in  the  following  sections.  The  Benchmark  Activity 
described  in  Section  3.2  focuses  primarily  on  model  verification  and  the  Precision  Test 
Modeling  program  described  in  Section  3.3  is  concerned  with  model  validation.  Both  of  these 
models  are  then  integrated  with  the  probabilistic  methods  described  in  Section  4  and  applied  to 
several  DNA  selected  problems  in  Section  5. 

3.2  BENCHMARK  ACTIVITY. 

The  goal  of  the  benchmark  activity  was  to  show  the  influence  of  calculational  approach  on  the 
predicted  response  of  a  circular  tunnel  in  a  jointed  rock  mass  subjected  to  ground  shock  loading. 
The  benchmark  problems  were  artificial  problems  that  tested  the  modeling  capabilities  and 
limitations  of  the  finite  element  and  distinct  element  computer  codes  used  by  the  participants. 
The  final  benchmark  calculation  was  preceded  by  several  preliminary  calculations  based  on 
simpler  but  related  problems.  The  common  element  in  the  calculations  was  the  geometric  and 
material  description  of  the  joints  and  intact  rock.  All  but  two  of  the  problems  solved  had 
analytical  solutions  to  compare  calculational  results  to. 

The  participants  in  the  Benchmark  Activity  were  RE/SPEC  (under  subcontract  to  SwRI),  R  &  D 
Associates  (RDA),  California  Research  and  Technology  (CRT,  and  now  TRT),  Weidlinger 
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Associates  (WA),  Itasca,  and  Lawrence  Livermore  National  Laboratory  (LLNL).  The 
benchmark  activity  has  been  well  documented,  so  only  a  brief  summary  will  be  presented  in 
here.  More  information  can  be  found  in  Simons  (1993),  Senseny  and  Simons  (1994),  and  Osnes 
(1991). 

To  summarize  the  effort,  several  bugs  and  modeling  errors  were  discovered  and  resolved  by  all 
of  the  participants  during  the  course  of  working  through  the  benchmark  problems.  The  results 
for  the  final  problem,  which  involved  a  deep  tunnel  in  an  explicitly-jointed  rock  mass  subjected 
to  stress  wave  loading,  three  of  the  participants  (CRT,  WA  and  Itasca)  obtained  (remarkably) 
similar  tunnel  deformation  results.  Using  PRONTO,  RE/SPEC  obtained  tunnel  closures  that 
were  substantially  higher  than  the  "consensus"  answer.  It  is  believed  that  the  larger 
deformations  were  due  to  the  explicit  joint  modeling  approach  described  in  Section  2.  However, 
because  an  analytical  solution  is  not  available  combined  with  the  fact  that  PRONTO  performed 
extremely  well  on  all  of  the  other  benchmark  problems,  some  question  remains  as  to  which 
answer  is  actually  correct. 

3.3  PRECISION  TEST  MODELING. 

3.3.1  Static  Tunnel  Tests  in  Limestone. 

In  many  model  validation  efforts,  predictive  models  are  simply  "calibrated"  to  agree  with 
experimental  results.  This  is  true  both  for  constitutive  models  calibrated  to  laboratory  data  and 
for  structural  response  models  calibrated  to  structure  test  data.  Because  of  this  practice,  the 
fidelity  of  the  model  is  unknown  when  used  to  predict  behavior  of  real-world  problems.  This 
necessarily  reduces  the  confidence  one  should  have  in  predictions  made  using  calibrated  models. 

One  way  to  assess  the  fidelity  of  a  predictive  model  is  to  apply  the  model  to  a  non-trivial 
laboratory  problem  that  is  representative  of  the  real-world  problem  without  any  prior  knowledge 
of  the  test  results,  and  then  critically  assess  the  performance  of  the  model  against  the 
experimental  results.  This  approach  is  the  one  being  taken  in  the  Precision  Test  Modeling 
(PTM)  program  sponsored  by  the  DNA.  The  goals  of  PTM  are  to  reveal  if  and  where  more 
research  and  development  is  needed,  explore  the  influence  of  computational  approaches  on  the 
outcome  of  the  numerical  simulations,  and  identify  where  modeling  approaches  differed  and  the 
relevance  of  those  differences.  The  motivation  for  such  an  effort  is  to  validate  the  ability  of 
current  state-of-the-art  predictive  models  to  simulate  real-world  problems.  The  PTM  program 
consists  of  laboratory  material  tests,  small-scale  structure  tests,  and  numerical  predictions.  Only 
the  numerical  predictions  are  summarized  here. 

The  PTM  participants  were:  Schwer  Engineering  &  Consulting  Services  (SE&CS  under 
subcontract  to  Aptek),  Southwest  Research  Institute  (SwRI),  Titan  Research  and  Technology 
(TRT),  and  Weidlinger  Associates  (WA).  The  static  and  dynamic  tunnel  tests  were  performed 
by  SRI  (Simons,  et.  al.,  1993). 

In  the  next  two  sections,  the  static  and  dynamic  tests  are  described  and  the  material  property  data 
used  in  the  exercise  tabulated.  Next,  the  calculational  methods  and  constitutive  models  used  in 
the  analysis  are  described  followed  by  comparisons  between  calculations  and  experiments. 
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Finally,  some  interpretation  from  the  comparisons  and  recommendations  for  further  research  are 
given. 

3.3,1 .1  Description  of  Static  Tests.  The  static  tests  consisted  of  cylindrical  test  specimens 
containing  lined  and  unlined  tunnels.  Several  load  paths  were  tested,  and  pressures,  strains  and 
tunnel  deformations  recorded.  From  the  suite  of  tests,  three  were  selected  for  analysis  in  the 
PTM  program,  herein  referred  to  as  STl,  ST2,  and  STS.  The  general  configuration  is  shown  in 
Figure  3-1.  The  test  specimens  were  all  30.5  cm  diameter  cylinders  taken  from  a  single  block  of 
Salem  limestone.  The  specimen  used  in  STl  was  30.5  cm  high  with  a  1.91  cm  unlined  tunnel 
spanning  the  diameter  of  the  specimen  at  midheight.  TTie  ST2  and  STS  tests  both  used  cylinders 
that  were  45.7  cm  high  (to  reduce  end  effects)  with  a  lined  tunnel  also  spanning  the  diameter  of 
the  specimen  at  midheight.  The  liners  were  fully  annealed  3003  aluminum  tubes  with  an  inside 
diameter  of  18.6  mm  and  a  wall  thickness  of  0.89  mm. 

The  STl  test  subjected  the  unlined  tunnel  to  a  single  uniaxial  strain  load-unload  cycle.  In  ST2, 
the  lined  tunnel  was  also  loaded  in  uniaxial  strain,  but  two  load-unload  cycles  were  performed. 
The  uniaxial  strain  path  is  of  interest  because  it  is  similar  to  the  loading  on  underground 
structures  from  detonations  above  ground  level.  In  ST3,  the  loading  approximated  a 
spherically-divergent  strain  path  similar  to  that  from  a  spherical  wave  propagated  by  an 
explosion  at  or  below  ground  level.  The  STl ,  ST2  and  ST3  loadings  are  shown  in  Figure  3-2 
(described  in  subsequent  sections). 


Figure  3-1.  SRII  static  test  configuration. 
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Figure  3-2.  Sandler-Rubin  cap  model.  Also  shown  are  the  load  paths  from  the  SRII  static  tunnel  tests  and 
the  smooth  cap  model. 


3.3.1 .2  Material  Properties.  A  significant  amount  of  test  data  for  Salem  limestone  has  been 
amassed  by  several  laboratories  working  independently  (ARA,  SRI  and  RE/SPEC  Inc.). 
RE/SPEC  was  tasked  with  performing  specific  tests  needed  to  characterize  the  parameters  in  the 
Sandler-Rubin  cap  model;  thus,  the  PTM  participants  were  directed  to  use  the  RE/SPEC  data 
(Possum,  et.  ai,  1995)  for  the  cap  model  parameters.  RE/SPEC  performed  a  large  number  of 
hydrostatic  compression,  triaxial  compression  and  uniaxial  strain  tests  on  Salem  limestone 
samples  taken  from  the  same  blocks  SRI  used  to  construct  the  static  and  dynamic  specimens. 
RE/SPEC  performed  two  different  fits  to  their  data.  The  first,  denoted  RE/SPEC-93,  was 
provided  to  the  PTM  participants  at  the  start  of  the  PTM  program.  The  second  fit,  denoted 
RE/SPEC-94,  included  lateral  strain  measurements  and  ultimate  shear  data,  presumably  to 
improve  the  fit.  The  RE/SPEC-93  data  was  used  for  this  study  because  the  RESPEC-94  data 
was  not  available  until  late  in  the  program.  For  completeness.  Table  3-1  lists  both  sets  of 
properties. 

The  lined  static  and  dynamic  tunnel  tests  used  aluminum  liners.  From  pressure  tests  performed 
on  the  liners,  SRI  measured  Young's  modulus  and  Poisson's  ratio  to  be  75.9  GPa  and  0.28 
respectively.  Since  significant  inelastic  deformation  on  the  liners  is  expected  in  the  tunnel  tests, 
SRI  also  measured  and  reported  the  effective  stress-effective  strain  response  of  the  aluminum. 
The  experimental  and  predicted  stress-strain  curves  are  shown  in  Figure  3-3. 
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Table  3-1.  Sandler-Rubin  cap  model  properties  for  Salem  limestone. 

Parameter 

Description 

RE/SPEC-93 

RE/SPEC-94 

Bulk  Modulus 

15725 

16214 

Shear  Modulus 

9112 

9531 

Failure  Surface  Parameter 

689.5 

199.43 

Failure  Surface  Parameter 

3.9e-4 

1.922e-3 

Failure  Surface  Parameter 

673.2 

185.77 

Hardening  Parameter 

1.44e-3 

8.616e-4 

w 

Hardening  Parameter 

0.08266 

0.108 

R 

Cap  Shape  Parameter 

4.215 

5.63 

(MPa) 

Cap  Location  Parameter 

-468.1  (-193.323) 

-441.26  (-131.217) 

Figure  3-3.  Comparison  of  experimental  and  predicted  aluminum  liner  response. 


3.3. 1.3  ralnilations.  Various  finite  element  discretizations  were  used  to  model  the  Static 

tunnel  tests.  SwRI  and  WA  used  full  three-dimensional  1/8  symmetry  models.  TRT  and 
SE&CS  used  "slab"  models,  which  neglected  the  outer  curvature  of  the  cylinder.  Since  the 
tunnel  responses  of  interest  were  at  the  center  of  the  cylinder,  the  error  in  neglecting  the  outer 
curvature  was  believed  to  be  small.  For  ST2  and  ST3,  ftictionless  contact  was  assumed 
between  the  liner  and  the  rock. 
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Responses  compared  included  axial  and  circumferential  strains  at  and  (see  Figure  3-1) 
and  springline  and  crown-invert  tunnel  closures. 

All  four  calculators  obtained  static  solutions  using  their  dynamic  codes:  SE&CS,  SwRI  and 
TRT  applied  the  loading  at  a  very  slow  rate  (quasi-static),  and  WA  used  dynamic  relaxation. 
The  calculations  were  initially  performed  using  a  two-invariant  cap  model,  and  comparisons 
between  calculators  made.  Since  all  except  SE&CS  used  exactly  the  same  two-invariant  model, 
this  allowed  a  consensus  check  to  be  made  before  employing  three-invariant  models.  The 
quasi-static  simulation  time  varied  by  calculator;  SwRI  used  20  ms. 

3.3. 1.4  Rftsiilts  and  Discussion.  Predicted  closures  from  the  STl  static  tunnel  simulations  are 
shown  in  Figures  3-4  and  3-5.  Using  the  two-invariant  models  (Figures  3-4a  and  3-4b), 
relatively  good  agreement  between  participants  is  seen  for  crown-invert  and  springline  closures. 
The  fact  that  SwRI,  TRT  and  WA  obtained  similar  results  is  because  they  used  the  same 
two-invariant  model.  SE&CS,  on  the  other  hand,  used  a  different  model  (see  Figure  3-2)  and 
obtained  different  results.  Comparisons  between  experiment  and  predictions  in  Figures  3-4a 
and  3-4b  show  considerable  disagreement.  The  crown-invert  closures  are  about  2%-3%  higher 
than  experimental  and  the  springline  closures  are  about  13%  lower  than  experimental. 

In  an  attempt  to  improve  the  closure  predictions,  three-invariant  yield  surface  models  were 
employed.  The  results  are  shown  in  Figures  3-4b  and  3-5b.  As  seen,  the  desired  effect  was 
obtained  —  the  springline  closures  increased.  However,  the  scatter  in  the  predictions  also 
increased.  Since  WA  and  TRT  results  compare  very  well,  and  also  since  they  both  used  the 
same  three-invariant  model,  the  scatter  between  all  participants  appears  to  be  due  to  the  use  of 
different  3-invariant  models.  This  indicates  that  closure,  or  more  generally  structure  behavior, 
is  very  sensitive  to  the  shape  of  the  yield  surface. 


a)  2-Invariant  b)  3-lnvarlant 


Figure  3-4.  STl  Crown-invert  closures:  a)  2-invariant,  b)  3-invariant. 
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Figure  3-5.  STl  Springline  closures:  a)  2-invariant,  b)  3-invariant. 


Results  for  the  ST2  and  ST3  tests,  both  of  which  include  liners,  are  shown  in  Figures  3-6  to  3-9. 
Figure  3-6a  shows  good  agreement  between  the  predictions  and  the  experimental  results.  The 
effect  of  including  the  third-invariant  was  small  as  shown  in  Figure  3-6b.  It  is  interesting  to 
note  that  SE&CS  predicted  very  little  recovery  during  unloading  (nearly  horizontal  lines), 
which  closely  resembles  the  test  data.  None  of  the  other  participants  predicted  this. 

Explanations  for  this  were  postulated  to  be  due  to  different  liner  constitutive  models  and/or 
closures  being  measured  on  the  rock  versus  the  liner.  Side  studies  performed  by  the 
participants  did  not  support  these  theories;  thus,  other  factors  may  be  to  blame,  such  as  the  use 
of  continuum  vs  shell  elements  to  model  the  liner. 

Springline  closures  for  ST2,  unlike  the  crown-invert,  showed  much  larger  variation  in 
prediction  when  the  third-invariant  was  included.  As  seen  from  Figure  3-7a,  predictions  made 
using  a  two-invariant  model  were  in  reasonable  agreement  with  experimental  results.  The 
inclusion  of  the  third-invariant  had  the  same  overall  effect  that  it  did  for  STl  —  springline 
closures  were  much  larger  and  the  scatter  between  predictions  increased.  Again,  this  scatter  is 
due  to  the  use  of  different  three-invariant  models. 
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a)  2-Invariant 


b)  3-Invariant 


50  100  150 

Vertical  Load  (MPa) 


200 


50  100  150 

Vertical  Load  (MPa) 


200 


Figure  3-6,  ST2  Crown-invert  closures:  a)  2-invariant,  b)  3-invariant. 


b)  3-Invariant 
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Figure  3-7.  ST2  Springline  closures:  a)  2-invariant,  b)  3-invariant. 


Crown-invert  and  springline  closures  for  STS  are  shown  in  Figures  3-8  and  3-9  respectively. 
Two-invariant  models  greatly  underpredict  both  the  crown-invert  and  springline  closures. 
Similar  to  STl  and  ST2,  the  crown-invert  closure  was  relatively  unaffected  while  the  springline 
closures  generally  increased  upon  inclusion  of  the  third  invariant.  SE&CS's  closure  results 
grew  very  large  because  their  calculation  predicted  that  stresses  in  the  free-field  (i.e.,  the  entire 
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specimen)  reached  the  failure  surface.  That  this  is  the  case  can  be  seen  in  Figure  3-2,  where  the 
STS  load  path  contacts  the  Smooth  Cap  (SM)  cap  failure  surface.  Their  closure  results, 
however,  should  be  simply  taken  as  qualitatively  "large"  since  complete  closure  and  specimen 
failure  will  occur  due  to  the  fact  that  the  failure  surface  is  perfectly  plastic.  The  fact  that  the 
STS  load  path  comes  very  near  to  the  Sandler-Rubin  (SR)  failure  surface  indicates  that  SwRI's, 
TRT's  and  WA's  closures  may  be  very  sensitive  to  small  changes  in  input  parameters. 


a)  2-Invariant 


b)  3-Invariant 


Figure  3-8.  ST3  Crown-invert  closures:  a)  2-invariant,  b)  3-invariant. 


a)  2-Invariant 


b)  3-Invariant 


Figure  3-9.  ST3  Springline  closures:  a)  2-invariant,  b)  3-invariant. 
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Also  worthy  of  note  is  the  fact  that  in  no  case  was  the  initial  slope  of  the  closure-load  curve 
predicted  well  by  the  models.  At  both  the  crown-invert  and  springline  location,  the  closures  are 
overpredicted  by  the  models.  Although  including  the  third-invariant  did  increase  the  final 
closure,  it  had  little  effect  on  the  initial  closures.  This  suggests  that  material  property  data  may 
be  incorrect,  or  other  mechanisms  not  included  in  the  models  may  be  involved,  such  as 
hardening. 

Two  observations  are  made  based  on  the  static  tunnel  simulations: 

1 .  Lx)ad-control  quasi-static  problems  can  be  extremely  sensitive  to  small  changes  in  the 
material  model  or  other  inputs.  Close  attention  has  to  be  paid  to  loading  rate  and  mesh 
discretization. 

2.  Static  problems  are  useful  for  verifying  dynamic  response  computer  models,  but  they  can  be 
highly-sensitive  to  small  changes  in  material  parameters,  and  can  introduce  additional  error 
if  true  quasi-static  conditions  are  not  achieved. 

3.3.2  Dynamic  Tunnel  Tests  in  Limestone. 

3.3.2. 1  Description  of  Dynamic  (SWAT-ID  Tests.  In  the  spherical  wave  tunnel  (SWAT)  test 
under  consideration,  a  spherical  explosive  charge  was  detonated  in  a  block  of  limestone 
containing  five  tunnels.  Recorded  data  from  the  experiment  included  particle  velocities  and 
stress  histories  at  several  locations,  and  post-test  tunnel  deformations.  The  general 
configuration  is  shown  in  Figure  3-10.  The  specimen  is  an  octagonal  block  of  limestone  1.2  x 

1.2  X  1.67  m  with  a  1 1 1-g  spherical  charge  of  pentaerythritol  tetranitrate  (PETN)  at  the  center 
of  the  block.  The  two  tunnels  selected  for  analysis  were  both  lined  and  located  at  distances  of 
14.5  cm  and  19.2  cm  from  the  center  of  the  PETN  charge.  The  tunnels  were  of  the  same 
diameter  and  liner  material  as  those  in  the  static  tests.  Stress  and  particle  velocity  gages  were 
located  as  indicated  in  Figure  3-10. 

Calibration  to  test  data  is  required  to  determine  either  the  modified  material  properties  used  by 
TRT  or  the  Duvaut-Lions  relaxation  parameter  x  used  by  SE&CS  and  SwRI.  Two  approaches 
were  taken:  1)  fitting  to  an  "average"  strength  increase  of  40%,  and  2)  fitting  to  the  SWAT-II 
free-field  velocity  or  stress  measurements.  The  40%  "average"  strength  increase  was  taken 
from  limited  testing  performed  by  ARA  on  Salem  limestone.  For  the  remainder  of  this  report, 
the  first  approach  will  be  referred  to  as  "uncalibrated"  and  the  second  approach  as  "calibrated." 
The  term  "calibrated"  is  used  because  a  portion  of  the  SWAT-II  experiment  (free-field 
response)  is  used  to  predict  the  tunnel  behavior,  i.e.,  the  predictions  using  "calibrated" 
viscoplastic  models  were  not  made  completely  independent  from  the  test  data.  The  values  for 
X  estimated  by  SE«&CS  and  SwRI  are  given  in  Table  3-2. 
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33.22  Calculations.  As  in  the  static  tests,  various  finite  element  meshes  were  used  to  model 
the  dynamic  tunnel  test.  SE&CS  and  WA  used  full  three-dimensional  models.  SwRI,  while 
originally  using  a  full  three-dimensional  model,  eventually  switched  to  a  "wedge"  mesh  to 
reduce  the  number  of  elements  and  speed  up  the  calculations,  shown  in  Figure  3-11.  TRT  also 
used  a  wedge-shaped  grid.  Frictionless  contact  was  assumed  between  the  liner  and  the  rock. 

Mesh  refinement  studies  were  performed  by  the  participants.  Consensus  reached  was  that  at 
least  five  rings  of  32  elements  around  the  circumference  of  the  tunnel  were  required  to 
adequately  model  the  complex  tunnel  wall  deformations  and  stress  gradients.  Beyond  this,  each 
participant  was  fi:ee  to  mesh  the  remainder  of  the  model. 


Figure  3-11.  SwRI  SWAT-II  wedge  mesh. 


In  Figure  3-10,  the  charge  cavity  radius  is  7.5  cm,  which  is  the  location  of  the  PV-1  velocity 
gage.  Thus,  the  detonation  of  the  PETN  was  not  simulated  in  the  calculations.  Instead  the 
recorded  velocity  transient  at  PV-1  is  used  as  the  driver  for  the  calculations. 

Responses  from  the  calculation  include  velocities  and  stresses  at  the  gage  locations,  and 
springline  and  crown-invert  tunnel  closures.  The  free-field  (no-tunnel)  environment  was 
computed  first  to  facilitate  comparison  between  calculators,  and  to  compare  with  the  free-field 
measurements  taken  in  the  experiment.  The  simulation  time  is  250/tf . 

3.3.2.3  Results  and  Discussion.  Predictions  of  free-field  velocities  were  computed  and 
compared  at  gage  locations  PV-2, 5, 8  and  11.  Using  rate- independent  models,  the  predictions 
agree  with  each  other  and  with  experimental  results  very  well.  The  only  exception  is  that  from 
PV-5  through  PV-1 1  the  predicted  peak  velocity  is  truncated.  However,  the  time-of-arrival, 
loading  slope,  positive  phase  and  negative  phase  all  agree  quite  well.  Employing  rate  models, 
SE&CS,  SwRI  and  TRT  were  able  to  improve  the  correlation  with  test  data.  The  same  general 
conclusions  can  be  made  with  stresses,  which  were  computed  and  compared  at  R-14.5  (see 
Figures  3-12  and  3-13)  and  19.2  cm.  In  general,  the  free-field  velocities  and  stresses  were 
found  to  be  much  less  sensitive  to  the  material  model  than  the  static  test  cases. 


Figure  3-12.  Free-field  radial  stresses  at  R=  14.5cm.  Predictions  based  on  rate-independent  model. 


a)  Uncalibrated  b)  Calibrated 


Figure  3-13.  Free-field  radial  stresses  at  R=  14.5cm.  Predictions  based  on  rate-dependent  models: 
a)  uncalibrated,  b)  calibrated. 
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While  rate-dependent  models  did  introduce  additional  scatter  in  the  predicted  results,  this  was 
outweighed  by  the  much  better  correlation  with  test  data.  As  will  be  discussed  in  the  following, 
this  was  not  generally  the  case  with  tunnel  closures. 

Tunnel  closures  for  the  R=14.5  cm  and  R=19.2  cm  tunnel  are  shown  in  Figures  3-14  to  3-17. 
For  the  R=14.5  cm  tunnel,  the  rate-independent  models  were  in  rough  agreement  (within  2%) 
with  experimental  data  at  the  crown-invert,  and  although  to  a  lesser  degree,  also  at  the 
springline.  Except  for  SwRI's  results,  which  agreed  with  experimental  results,  the 
rate-independent  models  underpredicted  crown-invert  closure  and  agreed  fairly  well  at  the 
springline  for  the  R=19.2  cm  tunnel. 

Both  SE&CS  and  SwRI  saw  very  little  change  in  their  crown-invert  closure  predictions  when 
rate-dependence  was  included  for  both  the  R-14.5  cm  and  19.2  cm  tunnels.  However, 
springline  closures  changed  dramatically.  The  inclusion  of  rate-dependence  changed  SE&CS' 
and  SwRI's  closures  from  2.5%  to  -1.5%  for  the  R-14.5  cm  tunnel  and  from  1%  to  -1%  for  the 
R-19.2  cm  tunnel.  These  predictions  actually  agree  with  the  experimental  (negative)  closures 
during  early  times  (i.e.,  springline  openture  is  measured).  Unlike  the  rate-dependent 
predictions,  however,  the  experimental  springline  opentures  quickly  reverse  and  become 
closure.  Using  their  rate-dependent  model,  TRT  predicted  higher  closures  at  both  the 
crown-invert  and  springline  locations  for  both  tunnels. 

SwRI's  and  TRTs  results  for  "uncalibrated"  and  "calibrated"  rate-dependent  models  are  also 
shown  in  Figures  3- 14b  to  3- 17b.  As  expected,  predictions  made  using  calibrated  rate  models 
were  closer  to  experimental  results;  however,  the  overall  quality  of  the  prediction  is  not 
improved  by  much. 


a)  Rate  Independent  b)  Rate  Dependent 


Figure  3-14.  SWAT-Il  crown-invert  closures  for  R-14.5  cm  tunnel;  a)  rate-independent,  b)  rate- 
dependent. 
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Crown-Invert  Closure  (%)  'g  Springline  Closure  (%) 


a)  Rate  Independent  b)  Rate  Dependent 
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re  3-15.  SWAT-II  springline  closures  for  R-14.5  cm  tunnel:  a)  rate-independent,  b)  rate-dependent. 


a)  Rate  Independent  b)  Rate  Dependent 


Figure  3-16.  SWAT-II  crown-invert  closures  for  R=19.2  cm  tunnel:  a)  rate-independent,  b)  rate- 
dependent. 
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a)  Rate  Independent 
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Figure  3-17.  SWAT-II  springline  closures  for  R-19.2  cm  tunnel:  a)  rate-independent,  b)  rate- 
dependent. 


By  including  rate-effects,  TRT  predicted  fundamentally  different  behavior  than  did  either 
SE&CS  or  SwRI.  TRT  used  independent  rate-enhancement  relations  for  the  failure  surface,  cap 
hardening,  and  three-invariant  yield  surface  relations.  SwRI  and  SE&CS  used  the 
Duvaut-Lions  linear  overstress  model.  Since  the  Duvaut-Lions  model  is  implemented  within 
the  constitutive  model,  it  accounts  for  the  spatial  and  temporal  variation  in  strain-rate.  In  the 
Duvaut-Lions  model,  the  yield  surface  is  effectively  moved  outward  (actually  the  viscid 
solution  lies  between  the  elastic  and  inviscid  cases),  which  explains  why  stresses  are  higher 
(less  yielding)  and  deformations  are  lower  (material  is  stronger).  Since  both  SE&CS  and  SwRI 
qualitatively  agree,  their  models  appear  to  be  consistent.  While  the  specific  details  of  how  TRT 
modified  their  cap  parameters  to  reflect  dynamic  effects  is  unknown,  this  approach  is  probably 
what  sets  TRT  apart  from  SE&CS  and  SwRI. 

3.3.3  Summary  of  PTM  Exercise. 

It  is  clear  from  this  exercise  that  more  research  is  needed  to  develop  "fully-defined"  material 
models  for  large-deformation  high  strain-rate  loading  of  geologic  materials.  For  instance, 
calculations  agreed  fairly  well  with  each  other  when  the  stress  states  were  predominately 
triaxial  compression,  but  differed  significantly  otherwise,  indicating  a  large  sensitivity  to  the 
shape  of  the  yield  surface  in  the  ii -plane.  Upon  further  investigation,  an  equally  large 
sensitivity  was  found  in  the  details  of  how  the  shear  failure  and  cap  (crush)  surface  intersection 
was  handled.  Large  variations  in  predictions  also  occurred  when  strain  rate  effects  were 
accounted  for. 

Comprehensive  testing  programs  should  be  carried  out  in  conjunction  with  numerical  modeling 
efforts  to  determine  which  effects  are  important  and  which  are  not.  For  example,  standard 
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material  tests  on  geologic  materials  primarily  probe  the  compressive  meridian  of  the  failure 
surface.  Through  comparisons  of  numerical  predictions  in  the  current  study,  it  was  found  that  a 
significant  portion  of  the  structural  response  was  due  to  yielding  closer  to  triaxial  extension 
than  compression;  however,  the  test  data  available  to  the  modelers  provided  little  guidance  for 
derining  the  model  in  this  regime.  The  next  logical  step  is  to  take  the  lessons  learned  from  the 
numerical  predictions  back  to  the  laboratory  and  obtain  more  test  data  and  repeat  the  modeling 
exercise.  Beyond  characterizing  the  location  and  shape  of  the  failure  surface,  test  data  is  also 
needed  to  characterize  joint  behavior,  material  damage,  and  flow  direction. 

Several  conclusions  from  the  PTM  study  can  be  made; 

1 .  The  constitutive  models  capture  the  main  effects  in  large-deformation  high  strain-rate 
loading  of  geomaterials;  however,  predictions  made  using  these  models  were  sensitive  to 
small  changes  in  the  model.  This  suggests  that  further  research  and  model  development  is 
needed. 

2.  More  comprehensive  and  possibly  non-standard  testing  is  needed  to  fully  define  the  material 
behavior.  The  testing  must  be  performed  in  conjunction  with  the  modeling  effort,  as 
opposed  to  being  performed  once  before  the  modeling  effort  begins. 

3.  Close  attention  must  be  paid  to  the  modeling  details  such  as  loading  rates,  grid  refinement 
and  boundary  conditions.  An  experienced  and  knowledgeable  analyst  is  perhaps  the  key 
factor  in  obtaining  high-quality  solutions  to  complex  highly-nonlinear  numerical 
simulations. 

The  objective  of  SwRI's  participation  in  the  verification  (Benchmark  Activity)  and  validation 
(Precision  Test  Modeling)  programs  was  to  develop  an  accurate  deterministic  numerical  model 
for  predicting  tunnel  deformation  under  high  stress  wave  loading  and  establish  confidence  with 
the  DNA  in  our  ability  to  apply  the  model.  The  approach  taken  was  to  enhance  the  PRONTO 
solid  dynamics  program  for  modeling  geomaterials  and  apply  it  to  realistic  problems  by 
participating  in  verification  and  validation  exercises.  The  results  indicate  that  solutions 
obtained  using  the  enhanced  PRONTO  model  are  quite  good  overall,  and  at  least  as  good  as 
most  other  numerical  models  currently  in  use.  This  provides  confidence  in  probabilistic 
predictions  made  using  the  enhanced  PRONTO  model. 
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SECTION  4 


PROBABILISTIC  ANALYSIS  METHODS 


4.1  INTRODUCTION. 

The  fundamental  problem  faced  in  probabilistic  structural  analysis  is  that  of  propagating  input 
parameter  uncertainties  through  a  numerical  model  to  compute  the  response  uncertainty.  Once 
the  response  uncertainty  is  known,  the  probability  of  failure  is  quantified  and  can  be  used  to 
make  decisions  regarding  design  changes,  etc.  Since  the  numerical  model  is  usually 
computationally-intensive  to  evaluate,  the  challenge  is  then  to  develop  a  probabilistic  analysis 
methodology  that  can  obtain  solutions  nearly  as  accurate  as  Monte  Carlo  simulation,  but  at  a 
much  more  reasonable  cost. 

Most  of  the  probabilistic  algorithms  described  in  this  Section  are  based  on  a  class  of  methods 
known  as  fast  probability  integration  (FPI),  which  have  been  shown  to  be  extremely  efficient 
relative  to  Monte  Carlo  simulation  with  comparable  accuracy.  For  probabilistic  finite  element 
analysis,  the  Advanced  Mean  Value  (AMV)  method  is  recommended.  Two  simulation 
methods — Latin  hypercube  and  Adaptive  Importance  Sampling — are  also  found  to  be  useful  for 
verifying  other  approximate  methods. 

4.2  UNCERTAINTY  MODELING. 

4.2,1  Random  Variables  and  Systematic  Errors. 

Input  uncertainties  arise  as  a  result  of  naturally  occurring  variations,  called  randomness,  in 
quantities  such  as  material  properties,  joint  spacing,  or  loading.  Another  form  of  uncertainty  can 
identified  as  modeling  or  systematic  error,  which  includes  errors  introduced  into  the 
calculations  by  the  use  of  simplified  models  (deterministic  and  probabilistic),  human  error,  and 
bias.  Randomness  is  intrinsic  to  nature  and  beyond  our  ability  to  control,  as  compared  to 
systematic  error  which  can  be  reduced  by  some  degree. 

Random  variables  are  modeled  using  a  joint  probability  density  function  (pdf).  The  joint  pdf 
provides  a  mathematical  relation  for  describing  how  the  input  parameters  vary.  Definition  of  a 
marginal  pdf  (i.e.,  for  a  single  variable)  is  typically  made  using  a  mean  value  /t,  a  standard 
deviation  o ,  which  is  a  measure  of  the  scatter,  and  a  distribution  type,  which  defines  how  the 
scatter  is  distributed.  As  an  illustration.  Figure  4-1  shows  the  pdf  for  a  variable  with  the  same 
mean  (30,000)  and  standard  deviation  (3000),  but  with  different  distributions  (normal  and 
extreme  value). 
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Figure  4-1.  Normal  and  extreme  value  pdf. 


To  determine  the  parameters  of  the  random  variables,  statistical  tools  such  as  the  graphical 
method,  method  of  moments,  method  of  maximum  likelihood,  least  squares  estimators,  and 
Bayes'  estimators  can  be  used  (Ang  and  Tang,  1984).  However,  in  most  cases  the  accuracy  of 
the  selection  is  not  clear,  especially  for  the  selection  of  the  distribution  type.  Therefore,  methods 
have  been  developed  to  test  the  quality  of  these  parameters  such  as  hypothesis  testing  and 
goodness  of  fit  tests. 

Approaches  for  modeling  random  variables  and  systematic  error  in  geomechanics  problems  are 
discussed  in  Tomg  and  Thacker  (1992b).  The  various  types  of  uncertainties  discussed  are 
broadly  classified  as  being  either  reducible  or  irreducible,  corresponding  to  systematic  error  and 
randomness  respectively.  The  approach  taken  is  to  use  a  probabilistic  approach  to  treat  both 
types  of  uncertainties.  Irreducible  uncertainties  are  modeled  as  random  variables,  which  results 
in  a  corresponding  uncertain  response.  Reducible  uncertainties,  on  the  other  hand,  are  treated  by 
assuming  Aat  the  mean  and  standard  deviation  of  the  reducible  random  variables  are  random, 
which  results  in  a  confidence  bound  on  the  calculated  probability.  Another  option  for  treating 
reducible  uncertainties  is  to  perform  two  (or  more)  probabilistic  calculations  to  bound  the  trae 
solution.  Although  the  form  of  the  result  from  bounding  analyses  bear  some  resemblance  to 
confidence  bounds,  the  manner  in  which  the  underlying  uncertainty  is  reflected  in  the  bounds  is 
different.  More  discussion  on  confidence  bounds  is  given  in  Section  4.8. 
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4.2.2  Random  Fields  and  Spatial  Variability. 

In  geomechanics,  spatial  variability  can  be  a  large  source  of  uncertainty.  Variables  possessing 
spatial  or  temporal  variability  in  addition  to  inherent  randomness  are  modeled  as  random  fields. 
An  example  of  a  random  field  is  the  spatial  variation  of  a  random  material  property  within  a 
layer  of  soil  or  rock.  Because  the  material  property  between  two  points  {i,j)  are  likely  to  have 
similar  values  if  they  are  located  close  to  each  other,  the  property  is  said  to  be  spatially 
correlated.  A  fully  correlated  random  field  implies  that  if  the  material  property  at  one  point  in 
the  layer  is  known,  all  other  points  in  the  layer  are  also  known.  At  the  other  extreme,  fully 
uncorrelated  random  field  implies  that  the  material  property  at  point  i  is  independent  of  point  j 
regardless  of  how  near  point  i  is  to  point  j .  In  most  cases,  the  random  field  is  partially 
correlated. 

Since  the  true  correlation  is  difficult  to  measure  between  every  point  i  and  j ,  correlation 
functions  are  usually  used  that  describes  the  manner  in  which  the  correlation  decays  with 
distance  or  time.  A  typical  correlation  function  is 

A 


where  A  is  the  distance  between  points  i  and  j  and  I  is  the  characteristic  length.  The 
characteristic  length  is  estimated,  measured,  or  related  to  other  characteristics  such  as  the  scale 
of  fluctuation.  Given  the  mean  and  standard  deviation  at  all  points  within  the  field  (usually  but 
not  necessarily  the  finite  element  nodes  for  spatial  fields)  and  the  correlation  matrix,  a  new  set  of 
uncorrelated  random  variables  can  be  obtained  by  transformation.  This  new  set  of  random 
variables  are  then  used  in  the  probabilistic  analysis.  The  transformation  is  exact  if  the  original 
variables  are  normally  distributed.  Further  details  on  modeling  random  fields  in  probabilistic 
finite  element  calculations  are  given  in  (SwRI,  1991). 

4.3  LIMIT  STATE  FUNCTION. 

Given  a  deterministic  tunnel  response  model  such  as  those  described  in  Section  2,  tunnel  response 
can  be  represented  in  the  form  of  a  response  function  such  as: 

Z(X)  =  Z(X,,X; . X,)  (4.2) 

where  Z  (X)  represents  the  tunnel  response  (e.g.,  closure)  and  X^  are  random  design  variables 
such  as  materiaJ  properties,  geometry,  or  loading.  To  indicate  failure,  a  limit  state  function  can 
be  written  in  terms  of  Eqn.  (4.2)  as 

g(X)  =  Z(X)-Z;,  =  0  (4.3) 

where  Zj,  is  a  value  of  Z(X)  such  that  g(X)  =  0  denotes  the  boundary  between  failure  and  non- 
failure.  Thus,  Eqn.  (4.3)  separates  the  design  space  into  "fail"  and  "safe"  regions.  The 
probability  of  failure  py  =  P  [g(X)  i  0]  is  defined  as  the  cdf  value  of  Z  at  Z^,  or 

p,=  (4.4) 
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where  /^(X)  is  the  joint  probability  density  function  for  X  and  Q  is  the  failure  region  given  by 
g  ^0.  The  complete  cumulative  distribution  function  (cdf)  (ranging  from  0%  to  100% 
probability)  can  be  constructed  by  simply  varying  Z„  and  re-evaluating  Eqn.  (4.4). 

Except  for  a  few  trivial  cases,  Eqn.  (4.4)  has  no  closed-form  solution.  Thus,  numerical  methods 
such  as  simulation  or  numerical  integration  must  be  used.  In  the  following  sections,  several  of 
these  methods  are  briefly  described. 

4.4  SIMULATION  METHODS. 

Simulation  methods  repeatedly  evaluate  the  deterministic  model  to  generate  a  sample  of  the 
response,  from  which  the  statistics  (e.g.,  mean,  standard  deviation,  etc.)  are  approximated.  A 
major  advantage  of  simulation  is  that  the  deterministic  model  does  not  have  to  be  simplified  or 
approximated  to  perform  the  analysis.  Also,  most  simulation  methods  are  exact  as  the  number 
of  simulations  approaches  infinity,  which  makes  it  attractive  for  checking  or  verifying  an 
approximate  method  of  probability  calculation.  The  main  disadvantage  to  simulation  is  that  the 
probabilistic  analysis  can  become  quite  costly  if  the  deterministic  model  is  complex.  Even  if  the 
deterministic  model  is  simple  to  evaluate,  the  analysis  can  get  time  consuminjg  if  small  failure 
probabilities  are  to  be  estimated,  which  is  usually  the  case.  This  can  be  demonstrated  using  the 
relation. 


%  error  =  200^  (4.5) 

y)  NPy 

where  pAs  the  failure  probability  and  N  is  the  sample  size  (Shooman,  1968).  There  is  a  95% 
chance  that  the  percent  error  in  the  estimated  probability  will  be  less  than  that  given  in 
Eqn.  (4.5).  Table  4-1  gives  the  total  analysis  time  for  various  probability  levels  and 
deterministic  function  evaluation  times  based  on  5%  error. 


Because  of  the  extremely  long  analysis  times  that  can  result  as  illustrated  in  Table  4-1,  many 
efficient  sampling  methods  have  been  developed.  Two  such  methods  are  described  in  the 
following  sections  after  a  brief  review  of  the  direct,  or  Monte  Carlo  method. 

4.4.1  Monte  Carlo  Simulation. 

Monte  Carlo  simulation  is  the  most  straightforward  simulation  method.  For  independent  random 
variables,  the  steps  involved  are: 

1 .  Generate  an  input  vector  by  sampling  from  the  input  variable  pdfs, 

X,  =  (4.6) 


where  ii  is  a  random  selection  from  a  uniform  distribution,  U[0,1],  and  FJ  ‘  is  the  inverse  cdf 
function  of  x. 

2.  Calculate  the  deterministic  response,  Z(X). 

3.  Evaluate  whether  or  not  failure  has  occurred,  g(X)  =  Z(X)-Z^  ^  0? 

4.  Repeat  steps  1-3  N  times, 

5.  Calculate  p^=  Ny/N,  where  Ny  is  the  number  of  samples  resulting  in  g(X) i  0. 

More  details  on  the  generation  of  independent  or  dependent  random  numbers,  mapping  from 
uniform  to  other  distributions,  and  the  calculation  of  response  statistics  are  given  in  (Ang  and 
Tang,  1984). 

4.4.2  Latin  Hypercube  Simulation. 

Monte  Carlo  simulation  can  yield  accurate  failure  probabilities  provided  the  number  of 
simulations  is  sufficiently  large.  Because  failure  probabilities  are  typically  small,  however,  the 
required  number  of  simulations  and  the  calculational  time  required  to  perform  a  single 
simulation  render  this  method  impractical  for  most  structural  reliability  applications  (see 
Table  4-1). 

Accuracies  commensurate  to  that  of  standard  Monte  Carlo  but  with  fewer  required  simulations 
can  be  achieved  with  constrained  sampling  schemes.  One  such  scheme  is  Latin  Hypercube 
simulation  (LHS)  (McKay,  et  al,  1979).  The  basic  idea  in  LHS  is  to  distribute  the  samples  such 
that  the  range  of  each  variable  is  efficiently  covered.  In  LHS,  each  of  the  k  random  variables  is 
divided  into  n  non-overlapping  intervals  based  on  equal  probability,  thereby  creating  k\n 
subregions  in  the  probability  space.  One  sample  is  taken  from  each  subregion  based  on  the 
probability  density  in  that  region.  Next,  the  n  values  for  the  first  random  variable  are  paired 
with  the  n  values  for  the  second  random  variable,  and  so  on  for  all  of  the  random  variables.  The 
n  random  vectors  of  length  k  are  then  used  to  perform  the  simulation. 
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The  pairing  of  the  variables  in  LHS  can  be  performed  in  such  a  manner  to  induce  a  user- 
specified  correlation  between  any  two  random  variables  (Iman  and  Conover,  1982).  The 
induced  correlations  are  based  on  the  non-parametric  technique  known  as  rank  correlation, 
which  remain  meaningful  in  the  presence  of  non-normal  distributions.  LHS  was  used  by 
Thacker  and  Senseny  (1994)  to  verify  a  probabilistic  cap  model  simulation  involving  correlated 
non-normal  variables.  This  will  be  discussed  further  in  Section  4.7. 

4.4.3  Adaptive  Importance  Sampling. 

The  approximate  response  function  around  the  MPP  (see  Section  4.5)  developed  using  the 
AMV+  (see  Section  4.6)  method  provides  only  approximate  probability  without  error  estimates. 
It  might  be  necessary  to  confirm  the  solution  with  minimal  extra  computations.  One  approach  is 
to  use  an  importance  sampling  method  that  performs  sampling  only  in  the  region  close  to  the 
MPP.  There  are  several  ways  to  develop  importance  sampling  regions.  In  the  adaptive 
approach  proposed  in  (Wu,  1994),  the  initial  sampling  region  is  generated  based  on  the  AMV+ 
function  and  gradually  increased  by  moving  the  sampling  boundary  until  the  sampling  region 
covers  the  failure  region  sufficiently. 

Figure  4-2  illustrates  a  curvature-based  adaptive  sampling  method,  in  which  the  sampling 
boundary  for  a  limit  state  is  a  parabolic  surface  that  is  rotationally  symmetric  about  the  vector 
OP  that  passes  through  the  origin  O  and  the  MPP,  P.  The  sampling  region  is  increased  by 
increasing  or  decreasing  the  curvature  of  the  parabolic  surface  until  there  are  no  more  failure 
points  in  the  increased  sampling  region.  Assuming  a  common  curvature,  the  samples  in  the 
sampling  region  are  generated  first.  Subsequently,  extra  samples  are  generated  only  in  the 
increased  region  introduced  by  the  curvature  change.  The  probability  increment  for  each 
curvature  change  can  be  evaluated  and  used  as  a  basis  for  generating  the  probability-consistent 
number  of  samples.  The  details  on  how  to  generate  samples  and  compute  the  failure  probability 
are  given  in  Wu  (1994). 


Figure  4-2.  Illustration  of  the  adaptive  importance  sampling  (AIS)  method. 
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The  generated  AIS  samples  can  further  be  used  to  compute  the  probability  sensitivities  with 
respect  to  the  probability  distribution  parameters  (i.e.,  mean  and  standard  deviation)  of  the  input 
random  variables.  The  computed  sensitivities  are  accurate  because  the  AlS-based  sensitivities 
are  based  on  the  "exact"  failure  regions.  Details  on  how  to  compute  probabilistic  sensitivities  are 
discussed  further  in  Section  4.9,  Wu  (1994),  and  Thacker  and  Wu  (1995). 

4.5  FAST  PROBABILITY  INTEGRATION. 

Z(X)  is  oftentimes  difficult  to  evaluate  requiring  a  computationally  intensive  methods  such  as 
finite  elements.  Thus,  the  key  to  an  efficient  probabilistic  analysis  is  one  that  minimizes  the 
number  of  times  Z(X)  must  be  evaluated.  As  an  alternative  to  Monte  Carlo  simulation, 
described  in  Section  4.4.1,  more  efficient  numerical  integration  methods  such  as  FPI  can  be  used 
(Wu  and  Wirsching,  1987;  Cruse,  et  ai,  1988;  Wu  and  Burnside,  1988).  The  entire  cdf  can  be 
efficiently  constructed  using  the  AMV+  (Section  4.6),  which  employs  the  FPI  algorithm.  A 
summary  of  the  basic  steps  involved  in  FPI  are  as  follows: 

1.  Obtain  an  approximation  to  the  exact  g-function  using,  for  example,  by  finite  differencing  or 
perturbation  methods. 

2.  Transform  the  original,  non-normal  random  variables  X  into  independent,  standardized 
normal  random  variables  (Ang  and  Tang,  1984),  u. 

3.  Calculate  the  minimum  distance,  p ,  from  the  origin  of  the  joint  pdf  (jpdf)  to  the  limit  state 
surface,  giO.  This  point,  u*,  on  the  limit  state  surface  is  the  most  probable  point  (MPP)  in 
the  transformed  u  space  because  it  is  the  most  probable  combination  (maximum  jpdf)  of  the 
standardized  normal  random  variables  on  the  g  ^  0  surface.  Figure  4-3  shows  a  jpdf  for  two 
random  variables,  the  limit-state,  and  the  MPP. 

4.  Approximate  the  g-function  g(u)  at  u*  or  g(X)  at  the  corresponding  X*  using  a  first-  or 
second-order  polynomial  function. 

5.  Calculate  py  =  P[giO]  using  available  analytical  solutions  (Madsen,  etal.,  1986;  Tvedt, 
1990). 
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Joint  Probability  Density  Function 
(JPDF) 


4.6  ADVANCED  MEAN  VALUE  METHOD. 

The  Advanced  Mean  Value  iteration  method  (AMV+)  was  developed  specifically  for  performing 
probabilistic  finite  element  analysis  (Wu,  et  al.,  1990).  The  algorithm  couples  the  finite  element 
code  with  an  FPI  code.  Two  different  AMV+  algorithms  have  been  developed:  One  alprithm 
is  used  when  the  value  of  the  cdf  is  prescribed  (i.e.,  when  the  probability  level  p  is  specified)  and 
the  corresponding  response  value  z  is  to  be  computed.  This  algorithm  is  termed  the  p-level 
procedure.  The  other  algorithm  is  used  when  the  z  value  is  prescribed  and  the  corresponding 
probability  is  to  be  computed  and  is  termed  the  z-level  procedure. 

The  p-level  procedure  is  used  to  compute  the  response  corresponding  to  a  specified  probability, 
or  p-level.  The  method  is  primarily  used  when  points  along  the  entire  range  of  the  cdf  are  to  be 
computed  (i.e.,  a  "cdf  analysis").  Table  4-2  gives  the  number  of  finite  element  solutions  that  are 
required  to  perform  a  first-order  probabilistic  analysis  using  AMV-i-.  In  Table  4-2,  N  is  the 
number  of  input  random  variables  and  M  is  the  number  of  p-levels.  MVFO  refers  to  "mean 
value  first  order"  and  AMVFO  refers  to  "advanced"  MVFO.  In  an  MVFO  analysis,  (N+1) 
structural  sensitivities  are  computed  and  used  to  construct  a  linear  performance  function  (g- 
function)  about  the  mean  values  of  the  random  variables. 

z®  =  ZW  +  t  HOC, 

=  Z,(X)  +  H(X) 
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In  step  2,  a  response  update  (or  "move")  is  performed,  which  entails  a  finite  element  solution  to 
update  the  response  value  at  each  of  the  M  probability  levels.  This  effectively  approximates 
H(X)  in  Eqn.  (4.7)  with  a  simplified  function  H(Zj).  If  the  change  in  Z  from  step  1  to  step  2  is 
larger  than  the  user-specified  tolerance,  Eqn.  (32)  is  re-computed  about  the  new  MPP  and  step  2 
repeated.  This  iteration  procedure  is  the  AMV+  method. 


Table  4-2.  Finite  element  solutions  required  to  compute  a  CDF  using  the  AMV+,  p-levels, 
procedure. 

Iteration _ SiSP _ Anqims _ Nhmbgr  pf  Finite 


0 

1 

MVFO 

N+  1 

0 

2 

AMVFO 

+  M 

(Move) 

1 

3 

First  Order 

+  (N+l)xM 

1 

4 

Move 

+  M 

2 

5 

First  Order 

+  (N+l)x  M 

2 

6 

Move 

+  M 

The  AMV+  z-level  procedure  is  used  to  compute  the  response  corresponding  to  a  specified 
probability,  or  z-level.  The  method  is  primarily  used  for  reliability  analysis,  where  p  must  be 
calculated  for  a  specific  value  of  z.  A  more  detailed  description  of  both  the  p-level  and  z-level 
AMV+  procedures  can  be  found  in  (Wu,  et  al.,  1990;  Thacker,  et  al.,  1991). 

4.7  CORRELATED  NON-NORMALS. 

If  the  original  random  variables  X  are  dependent,  a  transformation  must  be  made  to  a  space  in  which 
the  variables  are  independent  before  proceeding  with  the  FPI  calculations.  If  the  joint  pdf  is 
available,  the  Rosenblatt  (1952)  transformation  can  be  used: 

ttj  =  *-‘[F,(Xi)l 

“2  =  (4  8) 

!  ! 

=  O  UF,(xJx,.  Xj.  ....  X,.,)] 

where  represents  the  inverse  standard  normal  cdf  and  F,(*|‘)  denotes  the  conditional  cdf. 
Typically,  the  complete  joint  pdf  (jpdf)  is  very  difficult  to  obtain  based  on  either  theory  or 
experiment.  In  addition,  the  Rosenblatt  transformation  requires  the  inversion  of  the  conditional 
distribution  functions,  which  can  be  difficult  to  evaluate. 
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A  less  rigorous  but  more  practical  method  of  dealing  with  dependant  non-normal  variables  is  to 
approximate  the  jpdf  by  using  the  marginal  pdfs  and  associated  correlation  matrix 
(Hohenbichler,  et  al,  1981).  An  advantage  of  this  model  is  that  the  required  information  is 
relatively  easy  to  obtain  from  experimental  data.  The  first  step  in  the  procedure  is  to  transform 
the  dependent  non-normal  variables  X  to  standard  normal  u, 

u  =  <b-^[F^OO]  (4.9) 


Next,  the  correlation  coefficients  in  the  u-space  are  computed  using: 


/  / 


\  ^1  J 


^iui,updu,duj 


(4.10) 


where. 


= 


exp 


2(1  -  py 


(4.11) 


and  p„  „  is  the  unknown  to  be  determined.  This  model,  originally  proposed  by  Nataf  (1962), 
assumes' that  the  u-space  variables  are  jointly  normal.  As  a  result,  the  maximum  correlation 
assigned  in  the  X-space  can  be  limited  based  on  the  input  marginal  distributions  (e.g.,  two  non¬ 
normal  distributions  in  X-space  cannot  be  fully  correlated  and  also  be  jointly  normal).  This  is 
discussed  further  by  Liu  and  Der  Kiureghian  (1986),  who  also  present  approximate  solutions  to 
Eqn.  (4-10)  for  several  different  distribution  combinations. 

For  use  in  a  general  computational  framework  in  which  any  two  variables  may  be  correlated 
regardless  of  distribution  type,  a  more  general  solution  procedure  than  the  one  proposed  by  Liu 
and  Der  Kiureghian  (1986)  is  needed.  The  method  employed  here  approximates  Eqn.  (4-10) 
using  a  high-order  polynomial  expansion  and  the  normal  characteristic  function  (Wu,  et  al., 

1989;  Harren  1990). 

Next,  an  orthogonal  transformation  matrix  is  constructed  using  p  and  applied  to  transform  the 
correlated  normal  random  variables,  u ,  to  uncorrelated  normal  random  variables,  v, 

v  =  A^U  (4-12) 

In  Thacker  and  Senseny  (1994),  this  model  is  used  to  demonstrate  the  importance  of  including 
correlations  in  probabilistic  calculations.  A  limestone  test  specimen  is  subject  to  a  spherically- 
divergent  strain  path.  The  SR  cap  constitutive  model  is  used  to  model  the  material  response  m 
terms  of  nine  material  parameters.  Probability  distributions  for  each  parameter  and  the 
associated  correlation  matrix  were  measured  from  20  repeat  laboratory  tests  (Possum,  et  al., 
1995).  Probabilistic  calculations  were  performed  using  the  proposed  model  and  verified  using 
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Latin  hypercube  (correlated)  sampling.  The  results  show  the  approximate  method  described 
here  to  be  very  accurate  relative  to  Latin  hypercube  sampling.  The  results  also  show  significant 
differences  in  the  computed  probabilities  between  assuming  no  correlation  and  the  measured 
correlation,  thereby  indicating  the  importance  of  including  correlations  in  the  probabilistic 
calculations. 


4.8  CONnDENCE  BOUNDS. 

The  methods  for  probabilistic  analysis  described  in  the  previous  sections  all  assume  that  the 
parameters  of  the  random  variables  (mean,  standard  deviation,  etc.)  are  all  known  with  "high- 
confidence."  Statistically,  this  implies  that  a  large  amount  of  data  were  available  to  characterize 
the  random  variables.  In  most  applications,  however,  sufficient  data  is  not  available.  Thus, 
some  uncertainty  will  exist  in  the  parameters  of  the  random  variables  themselves. 

As  a  precursor  to  developing  an  efficient  procedure  for  constructing  confidence  bounds,  a  survey 
of  methods  for  uncertainty  modeling  in  geomechanics  is  given  in  Tomg  and  Thacker  (1992b). 
This  work  proposes  a  scheme  to  construct  confidence  bounds  for  structural  problems  containing 
uncertainties  due  to  inherent  variability,  estimation  error,  modeling  error,  and  human  error.  The 
effect  of  the  inherent  variability  is  reflected  in  the  calculated  reliability,  i.e.,  the  cdf,  whereas  the 
effect  of  estimation  error,  modeling  error,  and  human  error  are  reflected  in  a  confidence  bound 
on  the  calculated  reliability.  Strictly  speaking,  then,  the  bounds  should  not  be  termed 
"confidence  bounds,"  which  is  normally  used  to  denote  the  uncertainty  due  to  sampling  error 
only.  For  lack  of  a  better  name,  however,  "confidence  bounds"  will  be  used  here. 


The  proposed  nested  iteration  procedure  in  Tomg  and  Thacker  (1992a)  is  improved  by  Tomg 
and  Thacker  (1993)  to  utilize  die  probabilistic  sensitivity  factors  ( a )  computed  during  the 
AMV+  probabilistic  analysis.  The  basic  idea  is  to  constmct  a  first-order  Taylor  series  at  the 
point  of  interest  in  the  cdf 


ap 


(4.13) 


and  compute  the  derivatives  as  discussed  in  Section  4.9.  To  demonstrate  the  efficiency  of  the 
new  method,  the  procedure  is  applied  to  a  limestone  test  specimen  subject  to  divergent  strain 
path  loading.  Of  the  six  random  variables  considered  in  the  problem,  the  unconfined 
compressive  strength  and  joint  spacing  were  assumed  to  have  uncertain  mean  and  standard 
deviation.  Excellent  results  are  obtained  using  very  few  g-function  calculations. 

To  further  illustrate  the  method,  the  procedure  is  applied  to  the  analytical  tunnel  closure  model 
described  in  Section  2.  In  this  problem,  only  the  mean  of  the  external  pressure  is  assumed  to  be 
a  random  variable  with  its  own  mean  and  standard  deviation.  The  COV  for  the  mean  is  defined 
as  2%,  and  is  assumed  to  follow  a  normal  distribution.  None  of  the  other  random  variables’ 
parameters  are  modeled  as  random  variables. 

Figure  4-4  shows  the  cdf  and  confidence  bounds  for  the  50%  and  95%  confidence  levels 
computed  using  Monte  Carlo  and  the  method  described  here.  Any  confidence  level  can  be 
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chosen;  these  are  selected  only  to  demonstrate  the  type  of  result  obtained.  As  shown  in 

Figure  4-4,  the  proposed  method  obtains  excellent  results  at  a  fraction  of  the  cost  of  using  Monte 

Carlo. 


Figure  4-4.  50%  and  95%  confidence  bounds  for  the  analytical  tunnel  closure  model. 


4.9  SENSITIVITY  ANALYSIS. 

For  planning  and  design  purposes,  it  is  important  to  know  which  problem  parameters  are  the 
most  important  and  the  degree  to  which  they  control  the  design.  This  can  be  accomplished  by 
performing  a  sensitivity  analysis.  In  a  deterministic  analysis  where  each  problem  variable  is 
single-valued,  design  sensitivities  can  be  computed  that  quantify  the  change  in  some 
performance  measure  due  to  a  change  in  the  parameter  value. 

In  a  probabilistic  analysis  each  variable  is  usually  characterized  by  a  mean  value,  a  standard 
deviation,  and  a  distribution  type,  i.e.,  three  parameters  instead  of  one.  The  performance 
measure  is  usually  the  failure  probability  or  the  safety  index  p,  which  is  closely  related  to  the 
failure  probability.  Sensitivity  measures  can  be  defined  with  respect  to  each  of  these 
probabilistic  parameters  or  a  change  in  any  of  the  deterministic  parameters.  In  the  following 
sections,  the  calculation  of  deterministic  and  probabilistic  sensitivities  are  summarized  using 
both  sampling  and  numerical  probability  integration  methods. 
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4.9.1  Design  Sensitivities. 

Design  sensitivities  measure  the  change  in  response  due  to  a  change  in  an  input  parameter. 
Denoting  the  response  as  Z  and  the  input  variables  as  the  response  function  can  be  written  as 

Z  =  Z{X,,X^ . Z,)  =  Z(X)  (4.14) 


where  »  is  the  total  number  of  input  parameters.  From  this  the  i  -th  design  sensitivity  is  given 
by 

dZ 


4.9.2  Probabilistic  Sensitivities. 

Design  sensitivities  only  reflect  the  deterministic  importance  of  the  problem  variables.  In  a 
probabilistic  analysis,  meaningful  sensitivity  measures  are  how  py  changes  with  respect  to 
changes  in  input  parameters.  The  input  parameters  in  a  probabilistic  analysis  generally  include 
the  mean,  standard  deviation  and  distribution  type  for  each  problem  variable. 

The  limit-state  function 

g  =  Z(X)-Z  =  0  (4.16) 

is  written  such  that  for  a  specific  response  z,  {g:g  i  0}  denotes  the  failure  set.  Probabilistic 
analysis  methods  such  as  FPI,  FORM,  SORM,  etc.,  use  an  optimization  algorithm  to  locate  the 
minimum  distance  from  the  limit-state  to  the  origin,  denoted  p .  Calculations  are  typically 
performed  in  standard  normal  space,  or  b -space.  For  normally  distributed  variables,  the 
transformation  from  X  to  B  is  given  by  b,  =  (X^  -  /*,)/  o,.  For  non-normal  variables,  the 
transformation  is  usually  nonlinear  and  must  be  performed  numerically,  as  described  later. 

Once  the  minimum  distance  point  b/  is  located,  a  first-order  approximation  to  p^  can  be 
computed  from 

p,  =  ®(-P).  (4.17) 

where  *  is  the  standard  normal  cdf.  Thus,  sensitivity  measures  in  terms  of  p  can  also  be 
calculated  in  terms  of  p^  by  applying  Eqn.  4.17.  In  the  b -space,  P  is  given  by 

p  =  (4.18) 


so 


(4.19) 
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where  the  aj  are  the  probabilistic  sensitivity  factors.  FPI  methods  compute  a.\  as  a  by-product 
of  the  optimization  algorithm  from  the  g-function  gradient  calculation, 


«•  =  vg(«') 


N 


(4.20) 


which  show  that  the  a\  are  proportional  to  dgidui.  For  normally  distributed  variables, 
transforming  to  X-space  gives. 


iL  =  .g 

du;  dx;  dx; 


(4.21) 


which  shows  that  a]  are  proportional  to  the  design  sensitivity  dg/dX^  and  standard  deviation  of 
the  i  -th  variable.  For  non-normal  variables,  dg/dui  will  generally  also  be  a  function  of  the 
mean.  From  Eqn.  4.21  and  4.19,  it  is  seen  that  a\  reflects  the  i  -th  variable's  contribution  to  py. 

The  magnitude  of  a|  can  be  used  to  rank  the  random  variables  in  order  of  importance.  The  sign 
of  a]  can  be  used  to  determine  how  to  change  py  by  adjusting  aj .  Since  aj  is  a  function  of 
several  basic  parameters,  it  is  more  useful  to  know  how  py(or  P)  changes  with  respect  to  a 
change  in  the  mean  or  standard  deviation.  For  normally-distributed  variables,  these  are  easily 
derived  as 


li  = 

dp 

du; 

dui 

du; 

dui 

IP  = 

IP. 

du; 

do^ 

du; 

do, 

(4.22) 


For  non-normal  variables,  and  dtt,/do,  are  computed  numerically  as  follows.  First  the 

mapping  function  from  X-space  to  w-space  is  defined 

U,  =  <&-‘[F;(x^0)1  =  T(x,,e)  (4.23) 


where  is  the  cdf  of  x,  and  6  are  the  parameters  of  F^.  Next,  dUflddf  is  approximated  as 


du; 

di; 


t(x;.0)-t(x;,6) 


(4.24) 
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where  the  signifies  a  perturbation. 

So  the  procedure  is  to  perturb  6,  to  6^,  compute  T(x/,  6)  and  substitute  into  (4.24).  Although 
approximate,  only  the  mapping  is  changed,  so  no  further  g-function  calculations  are  required. 

4.9.3  Sampling  Based  Probabilistic  Sensitivities. 

The  samples  used  in  the  AIS  (or  Monte  Carlo)  analysis  can  be  used  to  compute  sensitivity 
measures.  Since  no  further  g-function  calls  are  required,  the  sampling-based  sensitivities  can  be 
computed  quite  economically.  The  basic  idea  is  presented  in  (Wu,  1994)  and  summarized  here. 
The  analysis  is  extended  to  include  lognormal  variables  and  applied  to  a  tunnel  closure  problem 
in  Thacker  and  Wu  (1995). 

By  definition 

^  "  L  "■  f  (4.25) 


where  Q  is  the  failure  region  and  is  the  pdf  of  x.  The  nondimensional  sensitivity  of  p  with 
respect  to  a  parameter  0  in  is  then 


dp/p 

ae/e 


r  •••  f 

/  &  =  E 

Iq  J 

480 

(4.26) 


where  E[*]  is  the  expectation  operator.  The  sensitivity  coefficients  S  with  respect  to  the  mean 
and  standard  deviation  are  thus 


dpip  _  g 

dpj  o, 

dp/P  _  E 

dOtlOi 

(4.27) 


In  (Wu,  1994),  5^^  and  are  derived  for  the  case  where  is  normally  distributed. 


nMonmU 

ftMormal 


(4.28) 
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The  sensitivities  can  also  be  computed  analytically  when  are  lognormal.  The  lognormal  pdf  is 
given  by 


/,  = 


E:q) 


1 

^  ln(x)  -  A,  ^  * 

'  2 

1  5  j  J 

(4.29) 


where  5 » \/ln(l  ♦  and  A,  =  In/i  1 1^/2.  Computing  df^l  and  df^^/dOf  and  substituting 
into  (4.27)  gives 


T” 


{«-«*- 1) 


2  ^ 


LV 


5 


ilogMimal  _  '^l 


uf-1 


-u. 


(4.30) 


where  =  (ln(ac)  -  A.)/? .  For  other  ^  /  dfif  and  ^  /  3o^  are  computed  numerically. 
4.9.4  Sensitivity  to  Distribution  Type, 


Another  important  sensitivity  is  how  py  changes  with  respect  to  the  distribution  type,  e.g., 
normal,  lognormal,  or  Weibull.  The  question  is  important  because  oftentimes  the  available 
experimental  data  is  adequate  to  characterize  the  mean  and  standard  deviation,  but  inadequate  to 
select  the  proper  distribution  type  with  a  commensurate  level  of  confidence.  The  choice  of 
distribution  type  becomes  increasingly  important  for  low  p^. 


The  sensitivity  measure  used  is  the  ratio  of  perturbed  ^  to  the  original  p 


(4.31) 


^  can  be  computed  by  re-running  the  probabilistic  analysis  with  the  distribution  type  changed 
for  one  of  the  random  variables,  or  estimated  using  an  approximate  procedure.  The  approximate 
procedure  involves  the  following  steps:  (i)  change  the  distribution  type  for  the  i-th  random 
variable,  (ii)  compute  a  perturbed  li/  using  the  new  distribution  type  and  (4.23),  (iii)  compute  a 
perturbed  ^  using  and  (4.18),  and  (iv)  compute  y  using  (4.31).  There  are  some  situations 
where  the  approximate  procedure  breaks  down.  For  example,  switching  from  normal  to 
lognormal  with  a  negative  variate  will  cause  the  procedure  to  fail.  However,  since  no  g-function 
evaluations  are  required,  the  approximate  procedure  can  be  an  effective  and  simple  method  for 
estimating  the  distribution  sensitivity. 
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SECTION  5 

NUMERICAL  EXAMPLES 

5.1  ANALYTICAL  TUNNEL  CLOSURE  MODEL. 

A  probabilistic  analysis  and  sensitivity  study  is  performed  using  the  analytical  tunnel  closure 
m^el.  This  model  describes  the  closure  of  a  circular  opening  in  a  Mohr-Coulomb  medium 
subject  to  uniform  internal  and  external  pressure  acting  on  the  opening  as  described  in  Section 
2.2.  The  problem  variables  are  listed  in  Table  5-1.  More  detail  on  the  probabilistic  analysis  can 
be  found  in  Thacker  and  Senseny  (1992). 


Table  5-1.  Random  variables  used  in  the  probabilistic  analysis. 

Variable 

Mean 

COV 

Distribution  Tvpe 

Young's  Modulus,  E 

30,000  MPa 

10% 

Extreme  Value 

Poisson's  Ratio,  v 

0.20 

10% 

Weibull 

Unconfined  Comp.  Strength, 

25  MPa 

20% 

Extreme  Value 

Friction  Angle,  <|> 

30° 

10% 

Normal 

Dilation  Coefficient,  x 

0.75 

5% 

Truncated  Normal 

0  X  ^  I'O 

External  Pressure, 

90  MPa 

10% 

Lognormal 

Internal  Pressure,  „ 

“a 

20  MPa 

5% 

Lognormal 

Using  the  AMV+  procedure,  the  cdf  (Figure  5-1)  and  a,  (Figure  5-2)  were  calculated.  As  an 
example,  for  vulnerability  assessments,  a  probability  of  kill  greater  than  0.85  is  desired.  This 

probability  corresponds  to  a  cdf  value  of  15%.  In  this  region  of  the  cdf,  two  variables  are 
dominate:  Young's  modulus  and  external  pressure. 
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Figure  5-1.  Cumulative  distribution  function  for  the  analytical  tunnel  closure  model. 


Random  Variable 

Figure  5-2.  Probabilistic  sensitivity  factors  (a)  at  cdf=15%  or  p^“0.85%.  As  shown.  Young  s  Modulus 
and  external  pressure  are  dominate. 
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Nondimensionalized  design  sensitivities  computed  as  part  of  the  AMV+  probabilistic  analysis 
are  shown  in  Figure  5-3.  E  and  Pj,  are  clearly  dominate;  however,  the  friction  angle  and  internal 
pressure  are  also  significant. 


E  Poisson  Sigu  Fricangl  Dilafact  Extpres  Intpres 
Random  Variable 


Figure  5-3.  Design  sensitivities  at  approximately  Pj^*»0.85. 


In  Figure  5-4,  the  probabilistic  sensitivity  with  respect  to  the  mean  (9p/9ft)  is  shown.  Since  P  is 
directly  related  to  probability.  Figure  5-4  indicates  that  a  change  in  the  mean  value  of  the 
dilation  factor  x  would  have  the  most  impact  on  the  probability. 


E  Poisson  Sigu  Fricangl  Dilafact  Extpres  Intpres 
Random  Variable 

Figure  5-4.  Probabilistic  sensitivity  with  respect  to  mean. 
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The  probabilistic  sensitivity  with  respect  to  the  standard  deviation  (dp/do)  is  shown  in 
Figure  5-5.  Here  it  is  seen  that  a  change  in  the  standard  deviation  of  Poisson's  ratio  would  have 
the  most  impact  on  the  probability. 


Figure  5-5.  Probabilistic  sensitivity  with  respect  to  standard  deviation. 


The  potential  error  in  assuming  the  wrong  distribution  type  for  the  random  variables  is  shown  in 
Figure  5-6.  The  sensitivity  measure  used  is  y  =  0/p  where  p  is  the  perturbed  beta 
corresponding  to  a  change  in  distribution  type.  Four  distributions  for  each  variable  are  tested. 
The  case  where  y  =  1  indicates  that  no  change  in  p  was  seen  when  the  distribution  type  was 
changed.  As  shown,  changing  the  distribution  of  E  from  it's  initial  value  (Extreme  Value)  to  any 
of  the  other  three  distributions  results  in  the  largest  change  in  p.  Poisson's  ratio,  friction  angle, 
external  pressure  and  internal  pressure  are  also  sensitive  to  changing  to  the  extreme  value 
distribution  in  particular.  In  an  actual  design  situation,  however,  certain  distributions  would  be 
ruled  out  based  on  physical  arguments  and/or  statistical  goodness-of-fit. 
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Figure  5-6.  Sensitivity  to  change  in  probability  distribution  type. 


Several  points  can  be  made  regarding  the  probabilistic  sensitivity  results.  First,  more 
information  is  made  available  for  decision  making  and  planning  than  is  available  from  a 
deterministic  analysis.  Only  certain  parameters  of  certain  variables  are  actually  controllable  so  it 
is  important  to  know  which  controllable  parameters  will  have  the  greatest  impact  on  the  design 
or  plan.  Second,  the  probabilistic  sensitivities  can  reflect  a  much  different  picture  than  design 
(deterministic)  sensitivities.  Because  of  this,  decisions  based  only  on  deterministic  analysis 
results  should  be  used  with  caution.  Finally,  from  a  computational  standpoint,  all  of  the 
sensitivity  results  derived  earlier  and  illustrated  here  are  very  economical  to  calculate;  most  are 
obtained  during  the  course  of  the  probabilistic  analysis.  For  those  that  must  be  computed  after 
the  probabilistic  analysis,  efficient  methods  as  described  earlier  have  been  developed  and 
implemented. 

5.2  PROBABILISTIC  DEEP  TUNNEL  ANALYSIS. 

The  probabilistic  response  of  a  deep  buried  tunnel  subjected  to  a  short-duration  stress  pulse  is 
demonstrated.  The  problem  is  highly  dynamic,  involves  nonlinear  material  behavior,  and  includes 
large  uncertainties  in  the  material  properties,  geometry,  and  loadings.  An  enhanced  version  of 
PR0NT02D  is  used  to  calculate  the  deterministic  tunnel  response.  The  failure  mode  considered 
is  excessive  closure.  The  AMV-i-  method  is  used  to  establish  the  cdf  of  tunnel  closure  with  relatively 
few  finite  element  re-solutions.  Monte  Carlo  and  AIS  methods  are  used  to  verify  the  probabilistic 
solution.  The  results  show  the  AMV-i-  method  to  be  highly  accurate  and  efficient. 
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The  problem  geometry  and  loading  are  shown  in  Figure  5-7.  The  analysis  considers  a  steel-lined 
tunnel  in  a  jointed  rock  mass  subjected  to  a  cylindrically-divergent  stress  wave.  Plane  strain 
conditions  are  assumed.  The  response  of  interest  is  maximum  tunnel  closure  at  any  point  in 
time. 


Pressure 


Figure  5-7.  Problem  description,  finite  element  mesh,  and  loading. 


The  elastic  jointed  rock  mass  is  modeled  using  the  compliant  joint  model  and  the  inelastic 
response  is  modeled  using  the  Drucker-Prager  model,  described  in  Section  2.3.3. 1 .  The  joint 
spacing  is  1  m  and  the  orientation  of  the  joint  sets  are  horizontal  and  vertical  with  each  set  being 
perpendicular  to  the  plane  of  the  problem.  The  steel  tunnel  liner  is  0.05  m  thick  and  is  modeled 
using  an  elastic-perfectly  plastic  material  model.  A  frictionless  contact  surface  is  assumed  at  the 
rock-liner  interface.  The  problem  variables  and  corresponding  values  are  listed  in  Table  5-2.  A 
preliminary  study  of  this  problem  using  fewer  random  variables  is  given  in  (Riha,  et  al.,  1992). 
The  rock  properties  are  typical  values  for  medium-strength  rock. 

The  deterministic  response  is  shown  in  Figure  5-8,  where  deformed  shape  plots  of  the  tunnel 
liner  and  adjacent  elements  for  several  times  through  the  analysis  are  shown.  All  tunnel  closures 
are  calculated  at  the  outer  surface  of  the  liner.  As  the  stress  wave  arrives  at  the  tunnel,  the  tunnel 
takes  on  an  oval  shape  with  closure  at  the  crown-invert  and  opening  at  the  springline.  Then,  as 
the  wave  passes  by  the  tunnel,  the  crown-invert  closure  decreases,  causing  the  springline  to 
begin  to  close.  This  is  evidenced  in  Figure  5-8  by  the  liner  pulling  away  from  the  rock  on  the 
sides  as  it  is  allowed  to  open  at  the  crown-invert. 
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Table  5-2.  Problem  variables  and  inputs. 


_ Variable _  Mean  Units _ COV _ Distribution 

Young's  Modulus _ 13.200  MPa _ 5^64 _ Extreme  Value 

Poisson's  Ratio _ 023 _ : _ 10 _ Weibull 

_ Peak  Load _ 100  MPa _ 10 _ Lognormal 

_ Density  (2) _ 2420  kg/m^3  0.2 _ Normal 

Unconfmed  Comn.  Strength  (21  51  MPa _ 06 _ Extreme  Value 


Friction  Angle  (1) _ 35  Degrees  14.3 _ Lognormal 

Dilatation  Angle  (1) _ 35  Degrees  14.3 _ Lognormal 

Tensile  Strength  2  MPa  0.00 _ : _ 


(1)  Fully  correlated.  (2)  Dropped  for  AMV+  analysis. 


Time  (ms) 


Figure  5-8.  Deformed  shape  plots  magnified  lOX  for  several  times  during  the  analysis  and  the  time 
history  of  crown-invert  closure. 
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To  save  computation  time  in  the  probabilistic  analysis,  the  deterministic  analysis  is  stopped  at 
150  ms,  when  the  loading  portion  of  the  stress  wave  had  passed  the  tunnel.  In  addition  to  an 
analysis  with  the  random  variables  at  their  mean  values,  two  perturbations  (±0.5a)  of  each 
random  variable  are  analyzed  for  a  total  of  15  deterministic  analyses.  The  additional 
perturbations  beyond  that  needed  to  construct  a  linear  response  function  are  performed  to  verify 
that  the  sensitivities  were  accurate  and  not  being  adversely  effected  by  the  unavoidable  noise 
associated  with  transient  direct  integration  solutions.  The  MV  cdf  based  on  these  15  FE 
solutions  is  shown  in  Figure  5-9. 

To  further  save  computation  time,  a  sensitivity  analysis  was  performed  based  on  the  MV 
solution.  From  this  analysis,  density  and  unconfined  compressive  strength  were  found  to 
contribute  only  a  small  amount  to  the  variation  in  crown-invert  closure.  Thus,  these  two 
variables  were  dropped  (i.e.,  held  at  their  deterministic  mean  value)  during  the  AMV-i-  analysis. 

To  calculate  converged  probabilistic  results  (MPP  locus),  the  AMV+  procedure  is  used.  Since 
the  primary  goal  of  the  analysis  is  to  verify  the  probabilistic  analysis  approach,  the  convergence 
tolerance  for  the  AMV+  iterations  is  set  to  1%,  which  is  smaller  than  what  is  typically  required. 
The  total  number  of  FE  solutions  to  compute  the  converged  cdf  using  AMV+  (including  the  MV 
solutions)  is  133.  After  confidence  is  gained,  the  convergence  tolerance  is  usually  relaxed  and 
only  single  perturbations  are  used  thereby  greatly  reducing  the  number  of  FE  calculations 
required. 

The  cdf  of  tunnel  closure  is  shown  in  Figure  5-9.  The  MV  cdf  is  accurate  at  the  50%  probability 
level  where  the  approximate  linear  response  function  was  initially  constructed.  However,  in  the 
tail  regions  the  difference  between  the  MV  and  AMV+  cdf  s  are  more  significant.  Even  though 
convergence  was  obtained  using  the  AMV+  procedure,  the  probability  calculation  is  still 
approximate  because  an  approximate  response  function  is  used.  To  check  quality  of  the 
probability  calculation,  Monte  Carlo  simulation  is  used.  As  shown  in  Figure  5-9,  the  results 
verify  the  AMV+  solution  from  about  -2a  to  +2o.  To  verify  the  solution  further  out  in  the  tails, 
the  AIS  method  is  used,  which  requires  on  the  order  of  100  FE  solutions  per  point  shown  in 
Figure  5-9.  To  verify  these  regions  of  the  cdf  with  Monte  Carlo  would  have  required  thousands 
of  FE  calculations.  It  should  be  noted  that  the  AIS  analysis  was  performed  only  to  verify  the 
AMV+  solution;  in  practice  a  larger  error  tolerance  for  AIS  would  be  used  resulting  in 
significantly  fewer  FE  solutions. 


56 


Figure  5-9.  Cumulative  distribution  function  for  crown-invert  tunnel  closure. 


The  probabilistic  sensitivity  factors  at  a  point  in  the  lower  and  upper  tails  are  shown  in 
Figure  5-10.  At  low  closure  levels,  Young's  modulus  contributes  the  most  to  the  closure 
variation,  whereas  at  higher  closure  levels  the  load  becomes  more  important. 


Figure  5-10.  Probabilistic  sensitivity  factors  for  the  deep  tunnel  problem. 
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This  analysis  demonstrates  that  advanced  probabilistic  algorithms  employing  fast  probability 
integration  methods  can  be  used  to  predict  the  uncertain  structural  response  of  underground 
structures  requiring  nonlinear  dynamic  calculations. 

5.3  TUNNEL  VULNERABILITY  ANALYSIS. 

The  problem  studied  is  a  shallow  tunnel  built  into  a  sloped,  multi-layered,  geology.  The 
conceptual  model  is  shown  in  Figure  5-11.  The  tunnel  is  of  size  2R  and  has  a  retaining  wall  of 
height  h  at  the  adit.  The  geology  considered  has  six  layers  of  varying  thickness  (constant  within 
each  thickness)  and  S-number  intersecting  the  tunnel  at  slope  0.  A  weapon  of  mass  hits  the 
hillside  at  velocity  V  at  a  distance  X  measured  up  from  the  retaining  wall  and  a  distance  Y  off 
the  axis  of  the  tunnel  and  penetrates  until  the  velocity  of  the  weapon  reaches  zero. 


Figure  5-11.  Conceptual  tunnel  vulnerability  model. 


The  following  assumptions  and/or  observations  are  made:  1)  the  layers  are  parallel  to  the  surface 
of  the  hillside  and  do  not  change  in  thickness,  2)  the  weapon  impacts  the  hillside  at  a  right  angle 
and  does  not  change  direction  while  penetrating,  3)  the  standoff  to  the  tunnel  is  the  absolute 
distance  from  the  tunnel  wall  to  the  center  of  gravity  of  the  weapon  after  the  weapon  reaches 
zero  (effective)  velocity,  4)  the  weapon  can  travel  through  the  tunnel  and  into  the  layers 
underneath  the  tunnel,  5)  and  the  weapon  may  not  intersect  the  tunnel  if  Y>R. 

A  table  of  the  model  inputs  is  given  Table  5-3.  The  mean  and  standard  deviation  for  s■^,  tj,  and  0 
were  provided  by  Chitty  (1995).  Data  for  the  remaining  problem  variables,  including 
probability  distribution  types,  were  assumed.  The  modeling  error  B  listed  in  Table  5-3  is  a 
multiplier  on  the  damage  frnction  in  the  probabilistic  calculation. 
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Table  5-3.  Model  variables  used  in  the  tunnel  targeting  problem. 


Random  Variable 


Layer  1  Thickness  _ 


Layer  1  S-number _ 


Layer  2  Thickness  _ 


Layer  2  S-number _ 


Layer  3  Thickness _ 


Layer  3  S-number _ 


Layer  4  Thickness _ 


Layer  4  S-number _ 


Layer  5  Thickness  _ 


Layer  5  S-number _ 


Layer  6  S-number _ 


Slope  Angle  _ 


Weapon  Initial  Velocity 


Transverse  Aimpoint  Error 


Modeling  Error 


Aimpoint  _ 


Retaining  Wall  Height 


Tunnel  Radius 


Weapon  Mass 


Damage  Function  Curve-fit  Coeffs. 


Identifier 


Mean  Value 

COV 

0.5  m 

25% 

10 

20% 

1.0  m 

25% 

6 

20% 

2.0  m 

25% 

2 

20% 

5.0  m 

25% 

1.5 

20% 

7.0  m 

25% 

1 

20% 

0.8 

20% 

33.81  degrees 

21% 

335  m/s 

10% 

0.0  m 

(s  - 1  m) 

1 

10% 

Variable 

0% 

3  m 

0% 

2  m 

0% 

8.9  kN 

0% 

Variable 

0% 

Probability 

Distribution 


Lognormal 


Lognormal 


Lognormal 


Lognormal 


Lognormal 


Lognormal 


Lognormal 


Lognormal 


Lognormal 


Lognormal 


Lognormal 


Lognormal 


Lognormal 


Normal 


Normal 


A  summary  of  the  deterministic  calculational  procedure  is  given  in  Figure  5-12.  For  a  given 
aimpoint,  the  first  step  is  to  compute  the  location  of  the  layer  interfaces.  The  tunnel  is 
considered  as  an  air  layer,  which  offers  no  resistance  to  the  weapon  should  it  be  encountered 
during  penetration.  As  shown  in  right  side  of  Figure  5-12,  the  presence  of  the  tunnel  could 
eliminate  a  layer,  or  the  tunnel  may  not  lie  on  the  penetration  path  for  some  values  of  ¥.  Next, 
the  penetration  depth  Zp  (measured  normal  from  the  surface  of  the  hillside  with  inward  being 
positive)  is  computed  using  the  PENCURV  penetration  program  (Adley,  et  a/.,  1994).  After  the 
penetration  depth  is  computed,  the  standoff  fi'om  the  weapon  center  of  gravity  (eg)  to  the  tunnel 
wall  is  computed.  Finally,  the  damage  (rubblized  volume)  is  computed  using  the  damage 
function,  which  is  assumed  to  be  a  curve-fit  to  several  known  data  points. 
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1 .  Calculate  interface  locations  (zi)  along 
line  of  penetration 

-zi  =  f(X,ti,  0,  h,  R,Y) 

-  Layers  extend  below  tunnel 

-  Layers  may  be  completely  masked  by 
tunnel 

-Tunnel  may  not  lie  on  penetration  path 

2.  Calculate  penetration  depth 

-Zp  =  PENCURV(zi,  Si,  Vp,  Mp) 

3.  Calculate  Standoff 

-  Sv  =  Abs[R  +  h  +  X  sin(e)  -  Zp  cos(0)l  -  R 

-  S  =  Sqrt[(Sv+R)^2  +  Y'^2)]  -  R 

4.  Calculate  Damage 

-  Vr  =  CO  +  Cl  (S)  +  C2(S)'^2  +  ...  +  Cm{S)''m 


(Weapon  also  tilted  out-of-plane) 


Figure  5-12.  Deterministic  calculation  procedure. 


Two  damage  functions  (Figure  5-13)  are  considered,  herein  referred  to  as  either  the  "long-tail" 
or  "short-tail"  damage  functions,  or/y  and)2  respectively.  The  short-tail  function  represents  an 
exact  fit  through  the  three  data  points  whereas  the  long-tail  function  represents  a  regression  fit 
through  the  three  data  points  and  another  point  at  {1 1 ,0}.  Both  functions  pass  through  {0,0}  and 
are  given  by: 


/l  - 


^  61.25  s  -  9.7  s'^  +  0.38 

0  <  S’  <11.454 

0 

othnwise 

34.76  s  +  1.67  -  0.71  s* 

0  <  s  <  8.24 

0 

otherwise 

(5.1) 
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Figure  5-13.  Damage  functions. 


BLU-1 09  into  multiple  soil  BH>1 09  into  multiple  soil 

layers  p<=0m)  layers  (X=1m) 


-1200b 


•laocb 


Y(mm) 


■ijnr' — ikiuj'  iW  9^ 

Y(mm) 


BLU-1 09  into  muittple  soil 
layers  (X=2m) 


BLU-1 09  into  multiple  soil 
layers  (X=3m) 


BLU-1 09  into  multiple  soil 
layers  {X=4m) 


BLU-1 09  into  multiple  soil 
layers  (X=5m) 


Figure  5-14.  Sequence  of  penetration  calculations  at  different  aimpoints,  X. 


As  an  illustration  of  the  deterministic  calculations.  Figure  5-14  shows  the  weapon  penetration  at 
six  different  aimpoints.  As  X  increases,  the  distance  from  the  surface  of  the  hillside  to  the  tunnel 
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increases.  Thus,  in  the  plots  in  Figure  5-14,  the  tunnel  gets  further  away  from  the  hillside 
surface  for  increasing  X.  Of  interest  is  that  for  X<3m,  the  weapon  passes  through  the  tunnel  and 
stops  shortly  after  reaching  the  other  side.  Since  the  eg  of  the  weapon  is  still  inside  the  tunnel, 
no  damage  is  done  as  modeled  by  the  functions  given  in  (5.1). 


/I 


0.7 

N 


g-6 

O 

Is 

c3 

®  4 

CL 


Standoff 


Penetration  Depth 


5  10  15  20  25 

Aimpoint,  X  (m) 


a) 


b) 


Figure  5-15.  Deterministic  calculational  results;  a)  variation  of  penetration  depth  (left  scale)  and 
standoff  (right  scale)  with  aimpoint;  b)  variation  in  damage  with  aimpoint. 


A  series  of  aimpoints  were  run  and  the  results  plotted  in  Figure  5-15.  Figure  5- 15a  shows  the 
variation  in  penetration  depth  and  standoff  with  aimpoint.  The  penetration  depth  increases 
initially  due  primarily  to  the  slope  of  the  hillside.  At  X®=2m,  the  penetration  depth  drops 
abruptly  due  to  the  fact  that  the  weapon  no  longer  penetrates  through  the  tunnel.  For  X>2m,  the 
penetration  depth  is  constant  and  the  standoff  increases  linearly  as  expected.  In  Figure  5- 15b, 
the  variation  in  damage  is  plotted.  For  both  either  the  long  and  short  tail  damage  functions,  the 
optimum  aimpoint  is  around  8  to  10m. 
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A  probabilistic  analysis  is  performed  using  the  deterministic  tunnel  vulnerability  model,  the 
inputs  listed  in  Table  5-3,  and  the  short-tail  damage  function.  The  cdf  of  damage  represents  the 
probability  that  the  damage  will  be  less  than  a  given  damage  value;  thus  the  vulnerability 
function  is  defined  as  1-cdf  such  that  vulnerability  denotes  the  same  meaning  as  probability  of 
kill  (pf).  For  example,  using  Figure  5-16  it  is  seen  that  pj^  is  approximately  90%  if  the  attack 
point  is  8m  and  the  kill  criteria  is  a  damage  level  of  100  m^.  In  other  words,  the  probability  that 
the  damage  will  be  greater  than  or  equal  to  100  m^  is  90%,  which  also  equals  p^.  can  always 
be  increased  by  decreasing  the  failure  criteria;  however,  as  shown  in  Figure  5-16,p;^  cannot  go 
above  90%  if  kill  is  assumed  to  occur  at  a  damage  level  of  100  m^.  From  Figure  5-16  it  is  also 
seen  that  goes  down  very  quickly  as  the  aimpoint  varies  from  8  m. 


Figure  5-16.  Vulnerability  function  (•sl-cdf). 


The  sensitivity  to  aimpoint  is  also  seen  in  Figure  5-17,  where  p/^  is  plotted  for  many  aimpoints. 
The  maximum  Pf^  is  found  (as  in  the  deterministic  calculation)  to  be  around  8  m,  but  the 
probabilistic  results  also  shows  that  pj^  drops  off  severely  as  the  aimpoint  varies  from  8  m.  The 
initial  peak  seen  in  Figure  5-17  corresponds  to  a  local  maximum  corresponding  to  the  weapon 
detonating  below  the  tunnel.  This  peak  was  not  found  in  the  deterministic  calculation.  To  see 
the  effect  of  choosing  different  kill  criteria,  two  kill  levels  in  addition  to  the  1(X)  m^  level  are 
plotted.  As  shown,  the  range  in  p^  is  about  20%  for  X=8  m  ±  5  m. 
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Figure  5-17.  Variation  in  probability  of  kill  with  aimpoint. 


Figure  5-18.  Nondimensionalized  sensitivity  of  with  respect  to  the  mean  of 
each  parameter. 
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Figure  5-19.  Nondimensionalized  sensitivity  of  with  respect  to  the 
standard  deviation  of  each  parameter. 


In  Figures  5-18  and  5-19,  the  nondimensionalized  sensitivity  coefficients  5^  and  5^  are  plotted 
(see  Section  4.9.3  for  the  definition).  measures  the  sensitivity  of  with  respect  to  the  mean 
value  of  each  random  variable,  and  measures  the  sensitivity  of  pj^  with  respect  to  the  standard 

deviation  of  each  random  variable.  Thus,  it  is  seen  that  changing  the  mean  value  of  the  S- 
number  for  layer  3  (S3)  and  thickness  of  layer  4  (t4),  or  the  standard  deviation  of  the  slope  angle 
will  have  the  largest  impact  on 

5.4  PROBABILISTIC  SWAT-B  ANALYSIS. 

The  refined  model  of  the  SWAT  n  dynamic  test  was  analyzed  probabilistically.  The  mesh  is 
shown  in  Figure  5-20.  The  deterministic  model  inputs  and  loading  are  the  same  as  those  given 
in  Section  3.3.2.  The  refined  model  used  here  was  created  after  the  verification  and  validation 
exercise  reported  in  Section  3.3.2  was  completed.  The  purpose  of  this  analysis  is  to  explore 
which  problem  variables  are  the  most  important  with  respect  to  tunnel  failure.  The  random 
variables  are  defined  in  Table  5-4  along  with  the  correlation  matrix  in  Table  5-5.  These 
probabilistic  inputs  were  obtained  by  Fossum,  et  ai,  (1995)  and  used  by  Thacker  and  Senseny 
(1994)  to  explore  the  significance  of  correlations  on  the  computed  cdf  results.  The  load  scale 
factor  term  listed  in  Table  5-4  acts  as  a  multiplier  on  the  velocity  time  history  used  as  input  to 
the  SWAT-II  calculation;  therefore,  it  is  given  a  mean  of  1.0.  The  standard  deviation  and 
distribution  were  selected  arbitrarily  to  be  0.1  and  normal  respectively.  There  is  no  correlation 
between  the  load  scale  factor  and  the  other  random  variables. 
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Figure  5-20.  Refined  mesh  used  for  the  probabilistic  SWAT-II  analysis. 


Table  5-4.  Parameter  summary. 


Load  scale  factor 


Mean  Value 

Std.  Deviation 

Prob.  Distribution 

15725 

761 

Lognormal 

9112 

290 

Lognormal 

689.5 

2.5 

Normal 

3.90E-04 

4.70E-06 

Lognormal 

673.2 

2.6 

Normal 

1.4EE-03 

8.49E-05 

Weibull 

0.08266 

0.00424 

Weibull 

4.215 

0.215 

Weibull 

-468.1 

4.5 

Normal 

1.0 

0.1 

Normal 

The  cdf  is  shown  in  Figure  5-21  along  with  the  50%  and  80%  confidence  intervals.  The  AMV+ 
analysis  procedure  is  used  to  compute  the  cdf  and  required  33  model  solutions.  The  confidence 
bounds  are  computed  based  on  the  coefficient  of  variation  (COV)  of  the  mean  and  the  standard 
deviation  using  20  test  samples  following  the  procedures  described  in  Section  4.8.  The  COV  of 
the  mean  and  standard  deviation  for  each  variable  is  shown  in  Table  5-6. 


1  Table  5-6.  COV  for  the  random  variable  mean  and  standard  deviation. 

Parameter 

Mean  COV 

Standard  Deviation  COV 

K  (MPa) 

0.019 

0.16 

G  (MPa) 

0.017 

0.16 

A  (MPa) 

0.019 

0.16 

B  (MPa'h 

0.428 

0.16 

C  (MPa) 

0.016 

0.16 

D  (MPa'h 

0.131 

0.16 

W 

0.131 

0.16 

R 

0.017 

0.16 

0.012 

0.16 

Load  scale  factor 

0.100 

6.10 
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Figure  5-21.  CDF  of  crown-invert  closure  showing  50%  and  80%  confidence  bounds. 


The  sensitivity  measures  are  the  change  in  Beta  with  respect  to  the  mean  and  standard  deviation. 
These  sensitivities  are  shown  in  Figures  5-22  and  5-23.  The  results  indicate  that  the  mean  value 
of  A  and  C  contribute  significantly  to  the  variations  in  tunnel  closure.  With  respect  to  the 
standard  deviation,  the  load  is  completely  dominate. 


Random  Variable 


Figure  5-22.  Nondimensionalized  sensitivity  of  |3  with  respect  to  the  mean. 


68 


Figure  5-23.  Nondimensionalized  sensitivity  of  p  with  respect  to  the  standard  deviation. 
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